# Using gdal_calc.py with conditionals on a float not giving correct results

I have two binary rasters (A and B) that I want to divide and get a binary raster (C) such that if A/B is greater than or equal to 50% (0.50), C = 1 and if not, C = 0.

Rasters A and B are of Byte type with unique values [0 1] (hence binary).

My command is as follows:

``````'gdal_calc.py -A rasterA.tif -B rasterB.tif --outfile=rasterC.tif --calc="((A/B) >= 0.50)" --NoDataValue=200 --type=Byte'
``````

The output raster looks good except that where pixels in rasterA = 1 and rasterB = 2, rasterC = 0 when it should be = 1 (1/2= 0.50). I'm guessing this is related to how in Python when you divide two integers by each other, the output truncates to an integer, but what do I need to change about the command to get my desired result?

Use gdal_translate to adjust your rasters to Float first and then you can do the math and avoid the integer division. http://www.gdal.org/gdal_calc.html

• Yeah but I thought only the intermediate calculation (A/B) needed to be float. The output will just be 1's and 0's and I want it to be byte. I'll try making the output of gdal_calc.py Float then doing gdal_translate to get it back to Byte. Thanks – user20408 Sep 19 '16 at 17:56
• Setting the output data type to Float still gives me 0.0000 where rasterA=1 and rasterB=2 .... Perhaps forcing either one of rasterA/rasterB or both be float would work. Similar to how we have to do float(A)/B or A/float(B) to get A/B as a float? I'm not sure – user20408 Sep 19 '16 at 18:08
• @user20408, I miss read your question, use gdal_translate to change the input rasters to float first then you can do the math. – Brian Sep 19 '16 at 18:11
• Thanks @Brian. Just writing raster A as a float did the trick. If you edit your answer I will accept it. – user20408 Sep 19 '16 at 18:36
• Great. Edited the response, sry about the misreading. – Brian Sep 19 '16 at 18:41

You can get what you want simply by using the following expression:

``````--calc="((1.0*A/B) >= 0.5)"
``````

There is no need to use gdal_translate to convert your source rasters as floats. Multiplying by 1.0 is casting the byte array as float, so the division is made in the float domain.

combine @Poiana Apuana and @Brian answers, I post below, hope its helpful. you just need to translate ONE of the image to float by `gdal_translate`, than use `gdal_calc.py` to get right results. Below is my full code to get BAIS2 Index：

``````import os
import sys
import glob
import numpy as np

def bais2(basedir):
"""
# param basedir
# Folder for storing S2 2 L2 products
"""
# Find out the images by tile
s2_folders = [item for item in os.listdir(basedir)
if not os.path.basename(item).endswith('BIAS')]

beta = 1e-4
# Traverse to find 20 meters res images with atmospheric correction done
for s2_folder in s2_folders:
B4_20m = glob.glob('{}\\{}\\**\\*B04_20m.jp2'.format(basedir, s2_folder), recursive=True)[0]
B4_20m_f = "{}\\BIAS\\floattmp\\{}_float.jp2".format(basedir, os.path.split(B4_20m)[1][:-4])
B6_20m = glob.glob('{}\\{}\\**\\*B06_20m.jp2'.format(basedir, s2_folder), recursive=True)[0]
B7_20m = glob.glob('{}\\{}\\**\\*B07_20m.jp2'.format(basedir, s2_folder), recursive=True)[0]
B8a_20m = glob.glob('{}\\{}\\**\\*B8A_20m.jp2'.format(basedir, s2_folder), recursive=True)[0]
B12_20m = glob.glob('{}\\{}\\**\\*B12_20m.jp2'.format(basedir, s2_folder), recursive=True)[0]
# B8_20m = glob.glob('{}\\{}\\**\\*B08_20m.jp2'.format(basedir, s2_folder), recursive=True)[0]

cmd1 = "gdal_translate -ot Float32 {} {}".format(B4_20m, B4_20m_f)
os.system(cmd1)
CMD = 'python P:\\AoiSoft\\ANACONDA\\Scripts\\gdal_calc.py' \
' -A {} -B {} -C {} -D {} -E {} --format=GTiff --outfile={}\\BIAS\\{}_bias2.tif ' \
'--calc="(1-1.0*sqrt((1.0*B*C*D)/(1.0*A*1e8+1e-9)))*((1.0*E-1.0*D)/(1e2*sqrt(1.0*E+1.0*D+1e-9))+1)"'\
.format(B4_20m_f, B6_20m, B7_20m, B8a_20m, B12_20m, basedir, s2_folder)
os.system(CMD)
pass

if __name__ == "__main__":
bais2("D:\\tmp")
``````
• Please, adjust all characters that are not English. – Taras Dec 5 '18 at 5:21