Help with rasterio.transfrom.xy


Cooney,Tom (ECCC)
 

Hi,

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[1], width=out_image.shape[2], count=1,

                                      dtype=rasterio.float64, transform=out_transform) as dataset:

                        dataset.write(out_image)

                        ds = dataset.read()

                        xs_to_mod = list(range(0, ds.shape[1]))

                        ys_to_mod = list(range(0, ds.shape[2]))

                        xs, ys = rasterio.transform.xy(out_transform, xs_to_mod, ys_to_mod, offset='center')

```
Thank you


Sean Gillies
 

Hi,

The xy method always returns projected coordinates. If you want long and lat, you must reproject the xy values using https://rasterio.readthedocs.io/en/latest/api/rasterio.warp.html?#rasterio.warp.transform.


On Tue, Jun 29, 2021 at 11:55 AM Cooney,Tom (ECCC) <Tom.Cooney@...> wrote:

Hi,

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[1], width=out_image.shape[2], count=1,

                                      dtype=rasterio.float64, transform=out_transform) as dataset:

                        dataset.write(out_image)

                        ds = dataset.read()

                        xs_to_mod = list(range(0, ds.shape[1]))

                        ys_to_mod = list(range(0, ds.shape[2]))

                        xs, ys = rasterio.transform.xy(out_transform, xs_to_mod, ys_to_mod, offset='center')

```
Thank you

,_

--
Sean Gillies