I try to take 1arc (1 arc-second res.) SRTM1 elevation data, project it into UTM, do slope, aspect and other calculations that require a meter-based projection, and finally reproject the results back to the original projection. At the moment I use a chain of gdal commands at the bash and a set of python files to achieve this, but is's messy and takes a lot of disk space.
So I plan to do all steps in a single script using rasterio. My basic idea was: - open LatLon source file with WarpedVRT, converting it to UTM on the fly - create multiple new raster bands (slope, aspect, other) that I add to the WarpedVRT instance - open a second WarpedVRT instance (this time projecting the extended WarpedVRT_UTM) to get my new data back to LatLon and the original resolution and extent
But by checking the array shapes I see the following: - original GeoTiff (LatLon): 3601x3601 px - WarpedVRT (UTM): 3904x3435 px [this is OK an expected] - WarpedVRT (LatLon): 3746x3788 px [this should be the original res and location]
I tried to get some inspiration from the commandline tools (basically window clip the 2nd WarpedVRT with something that resembles the '--like file' functionality of rio.clip but I'm stuck.
Any points if I am thinking along the right lines here? Or should I use another class?
My first attempts were along the following lines:
dst_crs = "EPSG:%d" % calc_utm_epsg(coord) with rasterio.open('s30_w072_1arc_v3.tif') as src: src_crs = src.get_crs() print src.shape # shape (3601, 3601) with WarpedVRT(src, dst_crs=dst_crs, resampling=Resampling.bilinear) as vrt_utm: print vrt_utm.shape # shape (3904, 3435) # (1) clone vrt for slope, aspect to hold results # (copy profile etc.), # maybe place new data bands into one vrt # # (2) do slope, aspect calulations # (...) # back-warp cloned vrt to original shape # change vrt_utm to slope_aspect_vrt with WarpedVRT(vrt_utm, dst_crs=src_crs, resampling=Resampling.bilinear) as vrt_src: print vrt_src.shape # shape (3746, 3788)
I'm using RasterIO v 1.0a9, GDAL2 on a Mac (homebrew)