[gdal-dev] reading/writing geometries using GDAL's MDA API

[Edzer Pebesma]

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 = _ ;
}