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)