Rasterio and WarpedVRT to go: Geo > UTM (compute stuff) > GEO without changing final extent?

Question

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)

Answer 1

If you create an instance of WarpedVRT without specifying height, width, and georeferencing, Rasterio calls a GDAL function to determine suitable values. It's quite possible that the function is not generally reversible. Instead, you should do this:

with WarpedVRT(vrt_utm, dst_crs=src.crs, dst_transform=src.transform, 
               dst_height=src.height, dst_width=src.width, 
               resampling=Resampling.bilinear) as vrt_src:
    # Do stuff.

That said, I'm not sure you'll be able to accomplish what you want to do with version 1.0 of Rasterio: the WarpedVRT object isn't backed by an XML file that you can edit and we don't have an API for adding bands to a WarpedVRT instance or modifying it in any other way. There is this trick for saving a WarpedVRT to disk, but that's it.