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): r = src.read(window=window) 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