1

I have a raster A .tif layer and a raster B .tif layer. I want to identify all pixels in raster B that are < 2 and set those corresponding pixels in raster A to 0, given that these rasters pixel grids are aligned. I understand that QGIS has conditional functions I can use for this within raster calculator, but since I need to loop through many raster files, I will need to rely and python and gdal. Trying to work this out, I would think code would look something like:

Calc("A where B < 2, 0, A", A=raster_A.tif, B=raster_B.tif, outfile=new_raster)

Though I am quite confused about using the proper logic and syntax here. How can I accomplish this conditional raster calculation using python and gdal?

1
  • are you bound to gdal or are you opento use other python libs (e.g. rasterio) ? Jul 18 at 7:23

1 Answer 1

1

I know that you are specifically requesting GDAL but I think rasterio could be easier to manipulate in your case and more pythonic. So if your environment authorize you to run pip install rasterio then here is a suggestion using pure rasterio + numpy trick:

import rasterio as rio
import numpy as np 

with rio.open("raster_A.tif") as ra, rio.open("raser_B.tif") as rb:
    
    # I assume the final raster will have the same characteristic as raster_A
    profile = ra.profile

    # load the data, I assume that you want to compare the 1st band of each file 
    raw_a = ra.read(1)
    raw_b = rb.read(1)

    # create the new dataset 
    data = raw_a * (raw_b >= 2)

    # make sure that data is the same type as the Tiff file
    data = data.astype(profile.dtype)

    with rio.open("new_raster.tif", "w" **profile) as dst:
        dst.write(data, 1)
6
  • Hello, rasterio is fine for me. Thanks! From what I understand in creating this new "data" raster is taking raster A, and from that recreating the same grid as all 0 values. And then on top of all those zeros, we take raster A and multiply it by all values from raster B that are >=2. And so for our "data" layer, if the raster B pixel is for example 3, that means the corresponding "data" pixel will be 3 * raster A pixel. And if that raster B value is 1.5, then the corresponding "data" layer pixel will be No_Data * raster A pixel? I am just a bit confused trying to understand the logic there! Jul 18 at 21:53
  • 1
    in numpy the result of a comparison is a boolean array so (raw_b >= 2) will be False (0) if raw_b is < 2 and True (1) if raw_b >= 2. as I'm multiplying it by raw_a, the final value of data will be raw_a where the condition is True, 0 elswhere. Jul 19 at 5:35
  • Ah ok, that makes a lot of sense now, thanks! Though in trying to run your code, on the line data = data.astype(profile.dtype) I received the error message: AttributeError: 'Profile' object has no attribute 'dtype'. I am a bit confused about what this line trying to do. Is it trying to determine if the raster is the same datatype as a numpy array? If "data" is a numpy array wouldn't it then not be the same data type as a tiff file? Might I be missing something there? Jul 19 at 20:50
  • 1
    I invite you to read the documentation of rasterio: rasterio.readthedocs.io/en/latest/# to better understand what each line is doing. the profile is embending the different characteristics of the Tiff file including the dtype but also the width, height, transform etc.... Jul 20 at 6:28
  • 1
    Ah I think I can see it now. This line seems to be setting the new raster to the same data type as raster A. And since "profile" here is actually a dictionary, I am accessing the "dtype" key to get the actual data type, which could be say "uint8" and then making sure the new raster is of that data type. At least that is what I can make of this. Thanks again for your help with this! Jul 20 at 22:38

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.