Is it possible to define a transform alongside a CRS, similar to geotiff?

[bretttully]

Let's say I have a geotiff with a well defined crs+transform. I then extract features from this raster.

I would like to save this data in the (pixel) coordinates that align directly to the geotiff. This has the nice property of keeping geometries specified in an integer coordinate system and means far lower likelihood of numerical issues in later processing that are otherwise introduced by round-tripping with floating point transforms.

e.g.

import geopandas as gpd
import rasterio
import rasterio.features
from shapely.geometry import shape

with rasterio.open("example.tif") as src:
    print(src.profile)
    layer = src.read(indexes=3)
    gdf = gpd.GeoDataFrame(
        [(v, shape(g)) for g, v in rasterio.features.shapes(layer)],
        columns=["value", "geometry"],
        # crs=???  # what should be set here?
    )
gdf.to_parquet("example-features.parquet")

The raster profile might contain something like

{
  'driver': 'GTiff',
  'dtype': 'uint8',
  'width': 2048,
  'height': 2048,
  'crs': CRS.from_epsg(4326),
  'transform': Affine(4.291534423828125e-05, 0.0, -75.5859375, 0.0, -3.296705832418409e-05, 39.84228602074339),
  ...
}

How would you put this into geoparquet? Is it even possible?

[brendan-ward]

Unless I'm missing something about what you are trying to do here, it should be this:

    gdf = gpd.GeoDataFrame(
        [(v, shape(g)) for g, v in rasterio.features.shapes(layer)],
        columns=["value", "geometry"],
        crs=src.crs
    )

The shapes extracted from the GeoTIFF are in the same CRS as the GeoTIFF.

From your transform, your pixels are not integer in size, so even if you converted your shape coordinates to pixel offsets within the GeoTIFF (translating them so that they are based on the origin point of the raster from the transform rather than the origin point of the CRS), they wouldn't be integer values, so you'd still have a floating point precision problem in addition to needing to define a custom CRS. Seems best to just keep your coordinates in the original CRS.

Feature coordinates don't have a concept of transform as used within a raster, which is really just a way of storing the parameters needed to calculate the position within the coordinate space specified by the CRS of each pixel within the raster.

[bretttully]

Perhaps I wasn't clear -- calling rasterio.features.shapes without specifying the transform kwarg means that the geometries are aligned to the pixels in the geotiff.

We can store these geometries into a regular parquet file with some metadata that describes the x, y, z, width, height of the original raster file. Downstream processing often doesn't care about the projection of the points (i.e. calculating intersections between features extracted from different layers), and so it is numerically nice if this is all performed in the same space of all values existing between 0, 2048.

However, it would also be nice if we could drag that parquet file into a tool like QGIS and have it land in the correct spot (like you can do with the geotiff, or a mapbox vector tile). For this, I guess we would need the transform just as the geotiff does.

Feature coordinates don't have a concept of transform as used within a raster, which is really just a way of storing the parameters needed to calculate the position within the coordinate space specified by the CRS of each pixel within the raster.

This probably answers my question... it is not something that is normally done.

[kylebarron]

You can store a per-row geotransform by storing an Arrow FixedSizeList with 6 float numbers per row.