Efficiently iterating through raster calculations with Rasterio or GDAL

I've got a very slow and round-about way to currently do this, but I struggle with the general concept of using windows in raster calculations and how to apply an affine or geotransformation between two large rasters. I'm hoping someone can help clarify what the heck i'm doing wrong here.

The following code is aimed at taking a large raster (200mb) and minusing a DTM from the raster values. The main points below:

I have these:

• DTM - NED raster tiles merged to bigtiff for the whole USA (77gigs)
• Water - Water depth raster over a portion of the US (water + DTM) (~200mb)

I want to perform a simple Water - DTM

Currently have a script that resamples a raster by individual columns where I query 1 point at a time to get DEM value at each point location in the raster. I know this probably isn't ideal but it's the only way I know to get a water depth raster efficiently.

The follow is my attempt at using Rasterio windows and then iterating through the pixels of each window to minus the DEM. At the end of the day i'd like to have another raster, but also a flat file with latitude/longitude/depth.

Most of this code came from here

with rio.open(tst) as src:
gt = src.get_transform()
a = src.affine
prj = Proj(src.crs)
a1 = a * Affine.translation(0.5, 0.5)
rc2en = lambda r, c: (c, r) * a1
nx, ny = src.width, src.height
print "Width: "+str(nx)+" | Height: "+str(ny)
for ji, window in src.block_windows(1):
cols, rows = np.meshgrid(np.arange(r.shape), np.arange(r.shape))
eastings, northings = np.vectorize(rc2en, otypes=[np.float, np.float])(rows, cols)
longs, lats = transform(prj, prj, eastings, northings)
for i, x in enumerate(list(r)):
if x > -100:
lat = lats[i]
lon = longs[i]
dem = xy2ras(lon, lat, dtmgt, bnd)
print "Lat: "+str(lat) + " | Long: "+str(lon)+" | Flood Height: "+str(x)+" | DEM: "+str(dem)
break
• If your Water raster and DTM have the same extents you can simply read each the data from each into a numpy array then subtract one array from the other. The resulting array can then be assigned to a new raster dataset. This eliminates iterating through every cell and will speed up your processing. – khafen Apr 9 '16 at 10:05
• So that's sort of my main issue here, I'm struggling with bounding the large DTM to the Water raster. I don't believe I can read the entire water raster into memory as an array, so what I've been trying to do is for each window of the water raster I read, I read only that bounding area from the DTM. – hydro_logic Apr 9 '16 at 19:52
• I think the first thing you should do is make sure both rasters are concurrent (see definition here etal.usu.edu/GCD/GCD5/GCD_GridConcurrency.pdf). This might mean you need to add empty cells to the water raster. Once you have that your rasters will share the same row and column indices, making it much more straightforward to perform the subtraction. Using python to iterate through each cell in a raster is slow, you'll need to set it up so you can read in each raster (or portions of each) as a numpy array to speed up processing – khafen Apr 11 '16 at 1:09
• I believe I have the concurrency issue under wraps right now. And I figured out that I could create a VRT of DTM to stack under the water raster. This means I can loop through windows and do a simple band arithmetic. Now I have to figure out how to write that back out, I'll update my question with the code when I finish it. Thanks for the advice! – hydro_logic Apr 11 '16 at 18:27