0

I have one GeoTIFF file with 2 bands. The data type is Byte. I want to calculate (Band2-Band1)/(Band2+Band1).

But the following command provided a dark image:

gdal_calc.py -A infile --A_band 1 -B infile, --B_band 0 --outfile outfile --calc "(A-B)/(A+B)" 

Does gdal_calc.py not calculate arithmetic with in the same file?

EDIT:

I run the following code using python:

import subprocess
gcalc = 'C:\\Users\\cleo\\Documents\\gdalpys\\gdal_calc.py'

a = 'D:\\a.tif'
b = 'D:\\b.tif'
outfile = 'D:\\result.tif'
expr = '(A-B)/(A+B)'
subprocess.call([sys.executable,gcalc,'--type','Float32','-A',a,'--A_band','1','-B',b,'--B_band','0','--outfile',outfile,\
    '--calc',expr],shell=True)  
a,b,outfile = None, None,None

gdalinfo -stats is as follows:

Driver: GTiff/GeoTIFF

Files: D:\tt.tif

Size is 36001, 36001

Coordinate System is:

GEOGCS["WGS 84",

    DATUM["WGS_1984",

        SPHEROID["WGS 84",6378137,298.257223563,

            AUTHORITY["EPSG","7030"]],

        AUTHORITY["EPSG","6326"]],

    PRIMEM["Greenwich",0],

    UNIT["degree",0.0174532925199433],

    AUTHORITY["EPSG","4326"]]

Origin = (-0.000138888888898,0.000138888888884)

Pixel Size = (0.000277777777778,-0.000277777777778)

Metadata:

  AREA_OR_POINT=Area

Image Structure Metadata:

  INTERLEAVE=BAND

Corner Coordinates:

Upper Left  (  -0.0001389,   0.0001389) (  0d 0' 0.50"W,  0d 0' 0.50"N)

Lower Left  (  -0.0001389, -10.0001389) (  0d 0' 0.50"W, 10d 0' 0.50"S)

Upper Right (  10.0001389,   0.0001389) ( 10d 0' 0.50"E,  0d 0' 0.50"N)

Lower Right (  10.0001389, -10.0001389) ( 10d 0' 0.50"E, 10d 0' 0.50"S)

Center      (   5.0000000,  -5.0000000) (  5d 0' 0.00"E,  5d 0' 0.00"S)

Band 1 Block=36001x1 Type=Float32, ColorInterp=Gray

  NoData Value=1.175494351e-038
  • try ruuning from the command line not from python (just to test) so you get the output and error messages printed to the terminal and edit your question to add them. – user2856 Jul 30 '14 at 2:17
  • yes I tested in terminal, and the error message i received is "the D:\\a.tif is not recognized....". Here is the input image I am working with commondatastorage.googleapis.com/earthenginepartners-hansen/… – cleo Jul 30 '14 at 2:51
  • You don't escape the backslash in the terminal Try D:\a.tif – user2856 Jul 30 '14 at 2:55
  • Sorry yes, actually I used only D:\a.tif in terminal – cleo Jul 30 '14 at 2:56
  • However I received the same error message 'D:\a.tif' is not recognized – cleo Jul 30 '14 at 3:16
1

This is an old question, but I've just run into what I think is a similar problem dividing two integer grids to get a percentage. The key for me was to convert them to float grids in the calculation. Just setting the output type didn't help. I am running GDAL 1.11.2. Example code is:

gdal_calc.py \
-A "_t2.tiff" \
-B "_t3.tiff" \
--calc="(A.astype(float)/B.astype(float))*100.0" \
--overwrite \
--type Float32 \
--outfile="_t4.tiff" 

A.astype(float) is numpy syntax to cast the integer grids to floats.

1

The output data data type defaults to the input data type which is Byte in your case. Your output needs to be floating point as all the values will be between -1 and 1. When the values get written to your output byte raster, they get truncated.

Use the --type option to set the output datatype.

gdal_calc.py --type=Float32 -A infile --A_band=1 -B infile, --B_band=0 --outfile=outfile --calc="(A-B)/(A+B)" 

The statistics of the output file should look something like this:

gdalinfo -stats outfile.tif
Driver: GTiff/GeoTIFF
Files: outfile.tif
       outfile.tif.aux.xml
<snip> coordinate system stuff </snip>
Band 1 Block=36001x1 Type=Float32, ColorInterp=Gray
  Min=-1.000 Max=1.000 
  Minimum=-1.000, Maximum=1.000, Mean=0.000, StdDev=0.480
  NoData Value=255
  Metadata:
    STATISTICS_MAXIMUM=1
    STATISTICS_MEAN=0.00044260699339504
    STATISTICS_MINIMUM=-1
    STATISTICS_STDDEV=0.47959178176737
  • Thanks for your response. Unfortunately, I am still getting the black image. – cleo Jul 29 '14 at 12:29
  • My full code is as follows: subprocess.call([sys.executable,gcalc,'--type','Float32','-A',infile,'--A_band','1','-B',infile,'--B_band','0','--outfile',outfile,'--calc','(A-B)/(A+B)'],shell=True) – cleo Jul 29 '14 at 12:33
  • @cleo Please edit your question and include the output of gdalinfo -stats for the output image. – user2856 Jul 29 '14 at 14:32
  • I have edited as per your requests – cleo Jul 30 '14 at 1:31
  • I am still facing this issue. All the solutions here didn't work for me. – wondim Jan 24 at 4:19

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.