4

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?

3 Answers 3

8

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.

2

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

5
  • 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, 2016 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, 2016 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, 2016 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, 2016 at 18:36
  • Great. Edited the response, sry about the misreading.
    – Brian
    Sep 19, 2016 at 18:41
0

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")
1
  • Please, adjust all characters that are not English.
    – Taras
    Dec 5, 2018 at 5:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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