Cropping and reprojecting a UTM raster to latlon

23 April 2024

GPXZ returns high-resolution rasters in UTM projection. This minimises reprojection on our end (as most of the GPXZ global elevation dataset is stored in UTM), and UTM projection lends itself nicely to analysis being in units of metres.

But some customers need results in a latitude & longitude coordinate reference system (mostly typically EPSG 4326), either for compatibility with other software, or to have the raster exactly cropped to bounds that are specified in lat/lon.

Take this bounding box

import rasterio

bounds = rasterio.coords.BoundingBox(-121.6271, 37.1622, -121.6184, 37.1650)

You can get a 1m resolution raster of the area using the GPXZ api:

resolution_metres = 1
url = f"https://api.gpxz.io/v1/elevation/hires-raster?key={GPXZ_API_KEY}&res_m={resolution_metres}"
url += f"&bbox_left={bounds.left}&bbox_right={bounds.right}&bbox_bottom={bounds.bottom}&bbox_top={bounds.top}"    

Rasterio can read these GPXZ urls directly, no need to save them to disk. We’ll use a WarpedVRT to convert the projection from UTM to latlon, then (despite only having a single raster) merge is the easiest way I’ve found in rasterio to crop.

import rasterio.merge

dst_crs = 'epsg:4326'

with rasterio.open(res) as f:
    f_vrt = rasterio.vrt.WarpedVRT(
        f,
        crs=dst_crs,
        resampling=rasterio.enums.Resampling.bilinear,
    )    
    zz, transform = rasterio.merge.merge(
        datasets=[f_vrt],
        bounds=bounds,
        resampling=rasterio.enums.Resampling.bilinear,
    )

and you’re done! zz is a numpy array of elevation data, in a lat/lon grid, aligned with the bounds initially passed to the GPXZ API.