I've been using GDAL's C++ multidimensional array (MDA) API lately to
read and write data cubes from R - excellent work! I was looking into
support for vector data cubes, multidimensional arrays with a single
dimension associated with a set of geometries.
What we can do is read and write dimensions, variables, and attributes,
but what (as far as I can tell) is missing is to read and write
OGRGeometry variables, using the CF specification for geometry [1],
which is supported by the NetCDF vector driver. Would it be feasible to
read and write CF-compliant geometries from multidimensional arrays
through the MDA API? Right now software writers have to sort this out
themselves, which is straightforward for points but cludgy for polygons.
Another issue is crs handling which I think could be done gracefully
using the OGRGeometry interface.
As a minimal data example: I added two ncdump files below for two time
instances and two POINT geometries. The first can be read and written by
the OGR API, but is a one-dimensional array that lacks the time
information of time series (time distributed over attributes/columns).
The second is a 2 x 2 multidimensional array with the time information.
GDAL's MDA API can read it but doesn't recognize geometries, GDAL's OGR
API cannot read it (obviously). I'd like to be able to read (and write)
the station dimension of the second one as an OGRGeometry, for points,
lines or polygons.
[1]
https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html#geometries
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
First example (written by OGR NetCDF driver):
netcdf t2 {
dimensions:
t2_node_coordinates = 2 ;
variables:
float t2 ;
t2:geometry_type = "point" ;
t2:node_coordinates = "t2_coordX t2_coordY" ;
double t2_coordX(t2_node_coordinates) ;
t2_coordX:axis = "X" ;
double t2_coordY(t2_node_coordinates) ;
t2_coordY:axis = "Y" ;
double t2_field_t1(t2_node_coordinates) ;
t2_field_t1:long_name = "Field t1" ;
t2_field_t1:geometry = "t2" ;
t2_field_t1:ogr_field_name = "t1" ;
t2_field_t1:ogr_field_type = "Real" ;
double t2_field_t2(t2_node_coordinates) ;
t2_field_t2:long_name = "Field t2" ;
t2_field_t2:geometry = "t2" ;
t2_field_t2:ogr_field_name = "t2" ;
t2_field_t2:ogr_field_type = "Real" ;
// global attributes:
:Conventions = "CF-1.8" ;
:GDAL = "GDAL 3.4.3, released 2022/04/22" ;
:history = "Tue Aug 16 12:14:25 2022: GDAL Create( t2.nc, ...
)" ;
data:
t2 = _ ;
t2_coordX = 0, 3 ;
t2_coordY = 1, 5 ;
t2_field_t1 = 0.33, 0.32 ;
t2_field_t2 = 0.04, 0.4 ;
}
second example:
etcdf t1 {
dimensions:
time = 2 ;
stations = 2 ;
variables:
double Precipitation(time, stations) ;
Precipitation:coordinates = "x y" ;
double stations(stations) ;
double time(time) ;
time:units = "days since 1970-01-01" ;
double x(stations) ;
x:axis = "X" ;
double y(stations) ;
y:axis = "Y" ;
double geometry ;
geometry:geometry_type = "point" ;
// global attributes:
:Conventions = "CF-1.6" ;
data:
Precipitation =
0.33, 0.04,
0.32, 0.4 ;
stations = 1, 2 ;
time = 19114, 19115 ;
x = 0, 3 ;
y = 1, 5 ;
geometry = _ ;
}
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
