Hi,
Not that I've a strong opinion on what GeoZarr/XArray should do, but yes
1D coordinate arrays are a bit of a pain to deal with, when they
actually just encode a regularly spaced grid. This is something I've hit
in the bridge between the GDAL multidimensional API (roughly netCDF
based) and the GDAL "classic" 2D raster API, when you want to expose a
view of a multi-dim dataset as a classic one. There's this bit of logic
in
https://github.com/OSGeo/gdal/blob/master/gcore/gdalmultidim.cpp#L7334
to guess if the 1D coordinate array is regularly spaced or not. The
netCDF driver has a similar heuristics too. This is quite error prone,
especially if Float32 coordinates are used. There it can be difficult to
detect if it is a regularly spaced array but Float32 hasn't sufficient
precision to fully encode it, or if it there are genuine irregularities
in the spacing. For Float64, the precision issues are less a practical
error, because comparing the expected value (assuming a regular spacing)
with the one of the 1D coordinate array with a very small relative
epsilon, like 1e-8, is sufficient.
Even
Le 03/04/2024 à 22:56, Michael Sumner via gdal-dev a écrit :
Here's an (ahem) extremely important discussion on the prospects for
xarray to extend from the coordinates-only model (like that of netcdf)
for geo reference:
https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/
I'm heartened to see this recognized at this level. I think GDAL per
se should feature in this discussion prominently, and not just the
downstream Python packages that are mentioned.
For my input, I'll be highlighting examples of
- degenerate rectilinear coordinates, and the problem of whether to
assign a regular grid there (with a trivial lightweight Translate a la
-a_ullr), or to define a new regular grid and push it through the Warp
api, as a dataset that picks up or can be pointed to the 1D coordinate
arrays. This has been the crux of the issue for me, I get to decisions
where it's not clear what was intended or what should be done going
forward.
- actual curvilinear coordinates, and the very very general power of
that to resolve to a regular grid with very simple specification
(specify some or nothing of target/extent/resolution/crs/dimension).
Key problems here are when the coordinate arrays are not auto-detected
(rare) or when the longitudes are unwrapped (less rare).
- that extent+dimension is complementary to the affine transform, and
way more user-friendly way to work with rasters for the assign and/or
target specification step (this is huge in R since the {raster}
package ca. 2009 and very well supported by Warp and Translate). It's
true that extent+dimension can't do shear params but it's a perfectly
valid way to work and flipping from transform to extent is trivial
numerically.
I hope some more voices from this community pitch in as well, very
happy to discuss here or offline about any of this. I'm not really a
technical expert but I have a good overview of the landscape that GDAL
sits in, and I'm making similar pitches for involvement in the R
communities.
Cheers, Mike
--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
