Cropping and reprojecting a UTM raster to latlon

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(url) 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.

The full code to then save that raster to disk would look like this:

import rasterio
import rasterion.merge


GPXZ_API_KEY = "your_api_key"

bounds = rasterio.coords.BoundingBox(-121.6271, 37.1622, -121.6184, 37.1650)
dst_crs = 'epsg:4326'
resolution_metres = 1
output_path = "latlon_raster.tif"

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}"  

with rasterio.open(url) 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,
    )

    profile = f.profile.copy()
    profile["crs"] = f_vrt.crs
    profile["height"] = zz.shape[0]
    profile["width"] = zz.shape[1]
    profile["transform"] = transform

with rasterio.open(output_path, "w", **profile) as f:
    f.write(zz)