plan for CNG 2025

[dcherian]

Hi Benoit

@maxrjones and I are both motivated to make some real progress on Xarray-for-geospatial in the next 3 weeks.

After a quick survey, I see at least two major items:

1. RasterIndex in rioxarray isn't moving. Can it be moved to this project ASAP instead? We can then move the generic affine pieces to Xarray later.
2. RasterIndex doesn't keep GeoTransform up-to-date after indexing. I'm happy to take on this update if you agree (xref this discussion)

With these done, a stretch goal would be to support more complex alignments (e.g. concat neighbouring tiles with same CRS and no overlap).

How does this sound to you?

[benbovy]

Hi Deepak,

@maxrjones and I are both motivated to make some real progress on Xarray-for-geospatial in the next 3 weeks.

Excellent!

RasterIndex in rioxarray isn't moving. Can it be moved to this project ASAP instead? We can then move the generic affine pieces to Xarray later.

I'd prefer limiting the scope of this project to (py)proj and CRS managing that is independent of the georeferenced data structure (raster, vector, irregular, etc.), if possible.

How about moving RasterIndex into a new repository? This would allow iterating quickly.

Regarding GeoTransform I see two options:

A. RasterIndex encapsulates an xproj.CRSIndex
B. RasterIndex only stores a (weak) reference to a spatial reference / grid mapping coordinate variable

I'd lean towards option B, which feels a bit hacky at first but has a few advantages:

• it keeps xproj.CRSIndex a 1st class Xarray index for any kind of geospatial dataset
• it will make the transition smoother if/when GeoTransform is decoupled from the spatial reference / grid mapping variable (hopefully), cf. Understanding concerns with CF encoding of CRS zarr-developers/geozarr-spec#20

We'd need to ensure that xproj.CRSIndex does not mess up with the GeoTransform attribute. It shouldn't (the default behavior of xarray.Index.create_variables() is to propagate the attributes of the input variables) but it'd require some testing. When the spatial reference variable is dropped, RasterIndex may simply ignore its (broken) reference and continue working as it does with its internal affine transform parameters.

What do you think?

[dcherian]

OK I completely agree with keeping the CRSIndex separate.

I'm not sure this project gets a lot more complex if you added RasterIndex but up to you...

For alternative projects, i can offer https://github.com/dcherian/rasterix/ . I am quite proud of the name, and it would be good if it were not just some random utilities :D

[scottyhq]

... crx was good, but rasterix is brilliant :)

[benbovy]

This is brilliant indeed!

I'll open a PR there with RasterIndex

I'm not sure this project gets a lot more complex if you added RasterIndex

Yeah indeed, but I thought Xproj as a lightweight Xarray wrapper around Pyproj would likely get adopted by various Xarray geospatial extensions (xvec, rioxarray, etc.)?

And Rasterix is definitely a good fit!

[dcherian]

I spent some time exploring this today (and then got sidetracked trying to get align to work hehe). See below for a list of workloads.

Basically I concluded that we need a layer that understands both the GeoTransform/Affine bit and the CRS bit. So I still lean towards (A). Though we could build an accessor/herlper layer that handled most of the complexity for the user.

Some "workloads" I thought through:

1. To go from spatial_ref with GeoTransform to a properly indexed dataset, we need quite a bit of setup.

   var = ds["spatial_ref"]
   x_dim = "x"
   y_dim = "y"
   
   from affine import Affine
   
   if geotransform := var.attrs.get("GeoTransform"):
      # TODO: delete GeoTransform if successful?
      geotransform = np.fromstring(geotransform, sep=" ").tolist()
   fwd = Affine.from_gdal(*geotransform[:6])
   centers = fwd * fwd.translation(0.5, 0.5)
   index = RasterIndex.from_transform(
      affine=centers,
      width=ds.sizes.get(x_dim),
      height=ds.sizes.get(y_dim),
      x_dim=x_dim,
      y_dim=y_dim,
   )
   
   newds = ds.set_xindex(["spatial_ref"], xproj.CRSIndex).assign_coords(
      xr.Coordinates.from_xindex(index)
   )

2. When available we should handle the GeoTransform with an AffineTransformIndex.
   a. We could probably delete GeoTransform from spatial_ref.attrs, and have the user manually write it back prior to writing to disk.
   b. It seems possible to create a GeoTIFF without GeoTransform but with CRS. Is that true? If so, what is the default Affine?
3. For xpublish-edr in particular, it'd be nice to have a .select(X=..., Y=..., crs=...) where X,Y are appropriate for the provided CRS. I can see this being useful even outside that project.

Something else we should think about is that odc-geo has this GeoBox concept that maps directly to RasterIndex.Affine. Perhaps this is another argument to keep it out of rioxarray, so it can be shared between the two subcommunities.

@benbovy Would 8.30AM MT, 4.30PM CEST on either Wed, Apr 9 or Thurs Apr 10 work for a coordination meeting? @scottyhq you're welcome to join too!

[dcherian]

To put this another way: what is "raster" about RasterIndex if it doesn't handle the CRS? At the moment, it's basically AffineTransformIndex.

[benbovy]

Though we could build an accessor/herlper layer that handled most of the complexity for the user.

Yes I think an accessor is the right layer for handling those things together.

Although I hope that with the coordinate transformation model proposed in zarr-developers/geozarr-spec#20 (comment) we could reduce workloads 1 and 2 to simply:

newds = (
    ds
    .set_xindex(["x", "y", "geo_transform"], rasterix.RasterIndex)
    .set_xindex("spatial_ref", xproj.CRSIndex)
)

Workload 2a seems fine to me, at least as a temporary solution. AFAIK other tools like rioxarray do not support automatic update of the GeoTransform attribute after indexing an object directly via Xarray's API. I thought it would be possible to have a constructor like RasterIndex.from_spatial_ref(variable) such that the index could continuously update the GeoTransform attribute of that input variable even if it is not associated with the index, but this will stop working as soon as the variable is copied. It is probably not worth it if the data model will be "fixed" in the mid-long term.

For workload 3 I was thinking that RasterIndex would basically be AffineTransformIndex + optional CRS (handled independently of the spatial reference variable) + xproj’s CRS-aware index interface (*). So workloads 1, 2 and 3 could be:

newds = (
    ds
    .set_xindex(["x", "y", "geo_transform"], rasterix.RasterIndex)
    .set_xindex("spatial_ref", xproj.CRSIndex)
    .proj.map_crs(spatial_ref=["x", "y", "geo_transform"])
)

newds.sel(x=…, y=…)  # CRS aware label-based selection

newds.proj("spatial_ref").crs == newds.proj("x").crs    # True

When there’s no spatial_ref or grid mapping variable it would be useful to allow passing the CRS directly as an option of RasterIndex:

newds = ds.set_xindex(["x", "y", "geo_transform"], rasterix.RasterIndex, crs=…)

---

(*) That’s actually the general model I see for Xarray geospatial indexes: data structure specific logic + optional CRS with custom API + xproj CRS-aware index interface. This model is also used by xvec.GeometryIndex (except currently for xproj’s interface, which would be easy to add). I can imagine other indexes be extended this way, e.g., for irregular grids or point clouds (xoak), staggered grids (xgcm), etc. Alternatively we could imagine each of those indexes optionally encapsulating an xproj.CRSIndex and thus be associated with a scalar grid mapping variable, but that would be less ideal IMHO.

---

It seems possible to create a GeoTIFF without GeoTransform but with CRS. Is that true? If so, what is the default Affine?

After reading the GeoTIFF specs it seems that the affine coordinate transform could be calculated from different combinations of GeoTIFF Tags (e.g., ModelTiepointTag for image corners + ModelPixelScaleTag instead of ModelTransformationTag).

AFAIK geo-referencing in GeoTIFF could also be done via a set of arbitrary ground control points + interpolation (the interpolation method not being part of GeoTIFF's specs). So maybe RasterIndex could eventually support that and may wrap something other than AffineTransformIndex?

For the default Affine, it seems that GDAL uses the identity matrix but it may not that simple from what I read in rasterio/rasterio#524.

[benbovy]

@benbovy Would 8.30AM MT, 4.30PM CEST on either Wed, Apr 9 or Thurs Apr 10 work for a coordination meeting? @scottyhq you're welcome to join too!

Either today or tomorrow 4.30PM to 5.00PM CEST would work for me.

[scottyhq]

Something else we should think about is that odc-geo has this GeoBox concept that maps directly to RasterIndex.Affine. Perhaps this is another argument to keep it out of rioxarray, so it can be shared between the two subcommunities.

Good point, and there are likely some differences to be aware of between rasterio/rioxarray and geo-box affine handling for certain cases (e.g. corteva/rioxarray#847)

It seems possible to create a GeoTIFF without GeoTransform but with CRS. Is that true? If so, what is the default Affine?

For the default Affine, it seems that GDAL uses the identity matrix but it may not that simple from what I read in rasterio/rasterio#524.

So if I understand correctly it can still be a 'valid' GeoTIFF to not have one but GDAL imposes (0,1,0,0,0,1)?

For workload 3 I was thinking that RasterIndex would basically be AffineTransformIndex + optional CRS (handled independently of the spatial reference variable) + xproj’s CRS-aware index interface (). () That’s actually the general model I see for Xarray geospatial indexes: data structure specific logic + optional CRS with custom API + xproj CRS-aware index interface.

I like the flexibility of this approach. I do find the setup cumbersome as well. Although users shouldn't have to do the setup - it can be handled by an open function/accessor right? (Basically I just want to do newds = rioxarray.open_rasterio('geo.tif') or rasterix.open('geo.tif').)

For what it's worth to my knowledge there are 4 primary ways rasters are georeferenced currently and an Affine is only one of them:

1. Affine
2. Ground Control Points (GCPs aka ModelTiepointTag)
3. Rational Polynomial Coefficients (RPCs)
4. External Geolocation Arrays (every pixel has a corresponding lon/lat, so basically CF model of explicit coords)!

Note that RPCs are not part of the GeoTiff Spec but RPCIndex could be really neat! In particular if you have them for an image, you can slice subsets and don't need to recompute the transform (see opengeospatial/geotiff#74). I've worked with GeoTiffs that actually might have all 3 approaches in the metadata (Affine, GCP, RPCs) and it's up to the user (via GDAL) which transform to use, but the default is to use the Affine if present.

[benbovy]

Although users shouldn't have to do the setup - it can be handled by an open function/accessor right? (Basically I just want to do newds = rioxarray.open_rasterio('geo.tif') or rasterix.open('geo.tif').)

Yes it should! Xarray IO backends (and/or eventually Xarray codecs) are a good place for doing this kind of setup using the metadata loaded from the file or store.

there are 4 primary ways rasters are georeferenced

!!

This would certainly justify maintaining raster index(es) and other utilities in a 3rd-party rasterix project! I also fully agree that this should be available to subcommunities like rioxarray, odc-geo, etc.

I think one major decision is "who handles the GeoTransform" attribute and keeps it up-to-date.

(@dcherian in an email conversation)

I agree, even more broadly I think the major decision is "who handles the spatial_ref / grid mapping variable". Should it be exclusively handled by xproj.CRSIndex? Or should it be bound to a geospatial index such as rasterix.RasterIndex, xvec.GeometryIndex, xoak.S2PointIndex, etc.?

I've been thinking more about it and TBH I'm still not sure I can come up with a detailed list of pros and cons for these two alternatives. Further collective thinking on this would be great.

Either way, I think it would be nice if we can stick with one of these alternatives across the whole Xarray geospatial ecosystem.

[dcherian]

Wow, thanks for the great discussion. There's a lot for me to learn and synthesize here, so lets meet tomorrow, April 10. I will send an invite.

[benbovy]

I agree, even more broadly I think the major decision is "who handles the spatial_ref / grid mapping variable". Should it be exclusively handled by xproj.CRSIndex? Or should it be bound to a geospatial index such as rasterix.RasterIndex, xvec.GeometryIndex, xoak.S2PointIndex, etc.?

Opened #25 for further discussing on this topic.

[benbovy]

Closing. The topics discussed here are being discussed elsewhere (in this repo or in rasterix).