3

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?

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

  • 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
5

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.

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

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.