We are working on a process that extracts values from raster data and returns them as GeoJSON geometries. Our use cases are point, line, and polygon.
For the line use case, we are trying to return the values and coordinates for each raster cell intersected by the line. We would also like to have the coordinates returned in geographic coordinates. To do this we are trying to use rasterio.transform.xy.
The problem we are encountering is that rasterio.transform.xy is not returning the coordinates in lat/long but rather in their native projection. We think this is caused by using the incorrect transform (Affine.affine) parameter.
Our setup is below. Note in this case shapes is a GeoJSON LineString:
with rasterio.open(raster_path) as src:
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
with MemoryFile() as memfile:
with memfile.open(driver="GTiff", height=out_image.shape, width=out_image.shape, count=1,
dtype=rasterio.float64, transform=out_transform) as dataset:
ds = dataset.read()
xs_to_mod = list(range(0, ds.shape))
ys_to_mod = list(range(0, ds.shape))
xs, ys = rasterio.transform.xy(out_transform, xs_to_mod, ys_to_mod, offset='center')