1

I am attempting to use the gdal_calc function to perform a calculation against 2 image bands.

I have a images a)B04.jp2 and b)B08.jp2. It is basically the red band and the NIR band. I want to perform the following: (b - a) / (b + a) Which is basically the formula to calculate the NDVI. I wrote the --calc parameter as follows: --calc="(B - A)/(B+A)"

So my entire call looks like this:

python .\gdal_calc.py -A C:\TMP\AWS\test\B04.jp2 -B C:\TMP\AWS\test\B08.jp2 --calc="(B - A)/(B+A)" --outfile='C:\TMP\AWS\test\test.tif'

Is this format correct? I received a blank image so I am thinking it may not be; however, I did receive a couple of

ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db errors while running the command.

I am pretty new to this and it took quite a bit to get this far.

Any thoughts?

Edit

Sorry, I should have guessed more information may be needed.

I installed Gdal through conda install via the conda-forge package. It is running version 2.3.3. I believe it is 64-bit based on the package name, py37hdf43c64_0.

I am running through powershell. I navigated to the Conda environment's script file, activated said environment, and ran the script exactly as pasted below.

Looking at it now, after a reboot, it appears as though some of the scripts are now missing. My gdal-info --formats is no longer working and I can't run it on the images either. It took me forever to get a conda install that had the jpeg2000 driver.

Any suggestions for getting a better, clean install on Windows bearing in mind my end goal would be to execute the gdal_calc script via command line?

EDIT 2

After a lot of struggling, I got the GDAL environment fixed, or at least I believe I did. Now I can provide the information requested.

I setup and installed everything to a virtual environment using Conda and the conda-forge channel packages.

python 3.8.2 h5fd99cc_4_cpython conda-forge gdal 3.0.4 py38h2fee047_1 conda-forge

Not sure based on the package information if it is 32-bit or 64.

I am now running the GDAL commands through the Conda command prompt, although I'd much prefer to do this via a standard PowerShell or Command Prompt window, as that would be the end goal.

Here is the GDAL info for raster A:

(GDAL) PS C:\> gdalinfo C:\TMP\AWS\test\B04.jp2 Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library Files: C:\TMP\AWS\test\B04.jp2 Size is 10980, 10980 Coordinate System is: PROJCRS["WGS 84 / UTM zone 16N", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["UTM zone 16N", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",-87, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",0.9996, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",500000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",0, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1]], ID["EPSG",32616]] Data axis to CRS axis mapping: 1,2 Origin = (300000.000000000000000,4500000.000000000000000) Pixel Size = (10.000000000000000,-10.000000000000000) Corner Coordinates: Upper Left ( 300000.000, 4500000.000) ( 89d21'52.85"W, 40d37'35.90"N) Lower Left ( 300000.000, 4390200.000) ( 89d19'50.34"W, 39d38'17.59"N) Upper Right ( 409800.000, 4500000.000) ( 88d 4' 0.53"W, 40d38'45.34"N) Lower Right ( 409800.000, 4390200.000) ( 88d 3' 5.22"W, 39d39'24.65"N) Center ( 354900.000, 4445100.000) ( 88d42'12.24"W, 40d 8'37.40"N) Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687 Overviews: arbitrary Image Structure Metadata: COMPRESSION=JPEG2000 NBITS=15

GDAL Info for raster B:

(GDAL) PS C:\> gdalinfo C:\TMP\AWS\test\B08.jp2 Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library Files: C:\TMP\AWS\test\B08.jp2 Size is 10980, 10980 Coordinate System is: PROJCRS["WGS 84 / UTM zone 16N", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["UTM zone 16N", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",-87, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",0.9996, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",500000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",0, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1]], ID["EPSG",32616]] Data axis to CRS axis mapping: 1,2 Origin = (300000.000000000000000,4500000.000000000000000) Pixel Size = (10.000000000000000,-10.000000000000000) Corner Coordinates: Upper Left ( 300000.000, 4500000.000) ( 89d21'52.85"W, 40d37'35.90"N) Lower Left ( 300000.000, 4390200.000) ( 89d19'50.34"W, 39d38'17.59"N) Upper Right ( 409800.000, 4500000.000) ( 88d 4' 0.53"W, 40d38'45.34"N) Lower Right ( 409800.000, 4390200.000) ( 88d 3' 5.22"W, 39d39'24.65"N) Center ( 354900.000, 4445100.000) ( 88d42'12.24"W, 40d 8'37.40"N) Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687 Overviews: arbitrary Image Structure Metadata: COMPRESSION=JPEG2000 NBITS=15

I am actually getting a black/gray/white image when processing now, but shouldn't I be getting a green/red/yellow?

And again, thanks a ton for the help, a lot of this is very foreign to me and I am trying to learn.

Improve this question
6
  • It sounds like a coordinate system or GDAL install problem, can you edit your question to include a GDALInfo of both rasters. Are you using 32 or 64 bit GDAL? What version(s) of python do you have? Did you install from the OSGeo or GISInternal installer (or are you not Windows)?
    – Michael Stimson
    Mar 9, 2020 at 23:52
  • Your GDAL_DATA and/or PROJ_LIB environment variables may not be set. Can you please edit your question and specify how you installed GDAL and exactly how you are running the python .\gdal_calc.py etc... (a standard cmd prompt, a cmd prompt started from a GDAL or OSGeo4W shortcut, a conda environment, something else...?)
    – user2856
    Mar 10, 2020 at 3:33
  • I think my environment became hosed after the system reboot last night. I have been having a terrible time attempting to get the gdal and rasterio suite loaded via Conda. Most of the time it will install but is missing scripts in the environment or it isn't installed with JPEG2000 support. Any suggestions for getting a clean, good install on Windows where the end goal would be to run the gdal_calc from command line?
    – Nathan Raley
    Mar 10, 2020 at 15:45
  • I just tried version gdal-300-1911-x64-core.msi from GS Internals and installing the Python 3.7 bindings GDAL-3.0.4.win-amd64-py3.7.msi for it, but there is no JPEG2000 driver listed. Any ideas?
    – Nathan Raley
    Mar 10, 2020 at 16:42
  • Make sure you tell it to use a float datatype if the input images aren't already: gdal.org/programs/gdal_calc.html#cmdoption-type
    – mikewatt
    Mar 10, 2020 at 19:54

1 Answer 1

Reset to default
2

About the installation of GDAL on Windows, use OSGeo4W.

About the grey output, NDVI pixels must have values between -1 and 1. If you want to render them with some color ramp you can use a Single Band Pseudocolor render style in a GIS application.

Improve this answer
2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.