[gdal-dev] discussion on xarray and the missing abstract-grid specification

[Even Rouault]

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: mdsumner at gmail.com
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240403/32e93b28/attachment-0001.htm>