3

I'm trying to calculate the "real height" of a region, meaning I want to sum the elevation value and building height value, if there is any building, for a region.

I know I can simply use Raster Calculator in ArcGIS and do a rasterA + rasterB, but I'm wondering if I can do it using Python.

I have generated two input raster files for this step:

First, Elevation raster

import rasterio

file1 = r"L:\lvpcgis\BowenLiuTemporary\ExtendLinesToLayer\Starry\allentown_elevation.tif"

elevation = rasterio.open(file1)

show(elevation, cmap = 'terrain')

It shows up like this:

Elevation

Second, Building raster with building height as pixel value:

file2 = r"L:\lvpcgis\BowenLiuTemporary\ExtendLinesToLayer\Starry\BuildingRasterTry2.tif"

bldg = rasterio.open(file2)

show(bldg)

It shows up like:

buildings

Now I'm wondering if I can add up these two layers so that the end result can reflect the sum of the values from the two rasters. It shall be doable but I haven't found any method online to achieve this with rasterio.

Other methods with Python will also suffice for me. I mentioned rasterio as it's the only module that I know can deal with raster files in Python.

Improve this question
0

Use raster.read(band_number) which returns a numpy array. numpy arrays can be summed together provided they're the same shape.

arr_bldg = bldg.read(1)
arr_elev = elevation.read(1)
arr_height = arr_elev + arr_bldg

https://rasterio.readthedocs.io/en/stable/quickstart.html#reading-raster-data

https://docs.scipy.org/doc/numpy/user/quickstart.html#basic-operations

Improve this answer
3
  • Thanks for the reply. My understanding is that rasterio will read rasters as a numpy array. But my two input datasets are apparently different in shape. Would it still work? – AndrewLebron Apr 25 '19 at 20:37
  • No, if they're different in shape numpy can't broadcast them because it's ambiguous which array elements correspond. Now you're getting into the need for windowed reads. Windows can be created from a bounding box expressed in map coordinates, so you can read the same geographic area from each raster. If I have time I can try to update my answer later – mikewatt Apr 25 '19 at 20:58
  • please sir if you can. I'm thinking if I can create an empty raster first with the size of the elevation raster and make all the value 0. And since the second raster is totally within the elevation raster, I can add these two to the empty raster and get the result this way. – AndrewLebron Apr 25 '19 at 23:47

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