1

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)

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.

  • Thanks for your comment. Bummer though that I cannot add bands to a WarpedVRT. Should I use a MemoryFile instead? Just guessing here... Any examples I could follow? 1) Load LatLon into MemoryFile 2) Reproject into another MemoryFile 3) Add bands from computation to 2) 4) Backproject 2) to 1) (same project/ same extent) 5) Write 4) to disk Is this the way to go? – Christian Sep 12 '17 at 19:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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