Unravelling a raster to a pandas dataframe

23 April 2024

Starting with this bounding box:

import rasterio

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

the GPXZ raster url will look like this:

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

After using rasterio to load the file, we’ll use the rasterio.Dataset attribute bounds and shape to define our pixel spacing. It’s important to use these attributes rather than the bounds originally used to query the elevation data, as GPXZ may buffer the output slightly, and our hires rasters are returned in UTM projection.

import numpy as np


with rasterio.open(res) as f:
    zz = f.read(1)
    x = np.linspace(f.bounds.left, f.bounds.right, f.shape[1])
    y = np.linspace(f.bounds.bottom, f.bounds.top, f.shape[0])

Finally, the numpy meshgrid function will compute coordinates for each pixel matching the elevation raster

import pandas as pd

xx, yy = np.meshgrid(x, y)
df = pd.DataFrame({
    'x': xx.flatten(),
    'y': yy.flatten(),
    'z': zz.flatten(),
})

Now you’ve got your dataframe, with one row for each point, with the elevation in column z. From here you can do other pixel-wise analysis, such as converting the coordinates to lat/lon:

from pyproj import Transformer

latlon_crs = 'epsg:4326'
transformer = Transformer.from_crs(f.crs, "epsg:4326", always_xy=False)
df['lat'], df['lon'] = transformer.transform(xx=df.x, yy=df.y)