1

I am looking for a way to perform raster calculations with two DEMs using rasterio. Basically, what the 'QGIS Raster Calculator' does. My problem is, that they do not have the same extent. Therefore, simple calculations with arrays don't work. Apart from that, they have the same CRS and resolution.

What i am looking for, is a way to expand the 'smaller' raster to the extent of the 'larger' so reading them into an array gives me overlapping arrays i can calculate with.

I've already looked into windowed reading, but that didn't get me were i wanted. Also building a VRT like suggested here could work, but i think there should be a more straight-forward way...

Two DEMs to perform raster calculations

Put very simply, here's what i want to do:

import rasterio as rio

dhm = rio.open(os.path.join(path_in, file_in_dhm))
tf = rio.open(os.path.join(path_in, file_in_tf))

arr_dhm = dhm.read(1)
arr_tf = tf.read(1)

result = arr_dhm - arr_tf

But that gives me - of course - an error:

operands could not be broadcast together with shapes (10561,16101) (7628,11309)

1
  • 1
    "I've already looked into windowed reading, but that didn't get me were i wanted"-- why not? That's how I'd do it
    – mikewatt
    Commented Nov 2, 2022 at 17:23

1 Answer 1

5

I think creating a window from the bounds of the larger raster then using that window with boundless=True to read the smaller raster should work.

Something like:

window = rasterio.windows.from_bounds(*larger_ds.bounds, transform=smaller_ds.transform)
smaller_array = smaller_ds.read(window=window, boundless=True) 
2
  • For some reason, when I do this and plot the updated 'smaller_array' equivalent, it looks like its extent has been properly expanded. However, when I go to add the two arrays together, it appears they are still not the same dimensions. Both raster datasets are the exact same projection. Not sure what is going on!
    – Ana
    Commented Aug 24, 2023 at 19:02
  • Set a fill_value perhaps? rasterio.readthedocs.io/en/latest/api/…
    – user2856
    Commented Aug 24, 2023 at 20:35

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.