Update GDAL NetCDF Vector Driver to support CF 1.8 geometries and Related Bugs / Feature Extensions

[wchen329]

The GDAL NetCDF vector driver [1] is currently compliant for the CF-1.6 convention. The CF-1.8 convention should be released soon with several new features including Geometry support [2], which NetCDF vector driver could greatly benefit from if made aware. I plan to update the NetCDF Read/Write Driver and possibly other files needed to support CF 1.8 geometries.

I'll update this issue and open a pull request when I have contributions ready for the GDAL community to review. Further requirements for such a pull request to be accepted (i.e. testing requirements) would be helpful.

[1] https://www.gdal.org/frmt_netcdf_vector.html
[2] https://github.com/cf-convention/cf-conventions/blob/master/ch07.adoc#geometries

Progress

• Baseline Support for Reading CF-1.8 geometries into OGRFeatures (Additions to netCDF Driver to allow Reading from CF-1.8 Encoded Geometries (#1287) #1541)
• Support for Writing from OGRFeatures into CF-1.8 compliant CF-netCDF Files (netCDF CF-1.8 Driver Updates - Writing Support, Bug Fix, Streaming etc. #1680)
• Support for CF-1.8 Style Grid Mappings and Geo-References (netCDF CF-1.8 Driver Updates - Writing Support, Bug Fix, Streaming etc. #1680)
• Bug: NetCDF-4 files (FORMAT=NC4) work for writing CF-1.8 point geometries, but otherwise fail for the other feature types. Resolve by: using infinite dimensions instead of fixed dimension size for NC4 specifically and making other changes
• Bug: Update GDAL NetCDF Vector Driver to support CF 1.8 geometries and Related Bugs / Feature Extensions #1287 (comment) (Yahara CRS and something with States). [PR]
• Feature Request: Ability to remove dimensions M (in general with warning) or Z (if not particularly lossy?). [PR]
• Optimization: Do not pre-define dimensions for NC3, only define them upon first write.

[rouault]

CC @tomkralidis Perhaps you have some comments on this ? I'm not sure how critical it is to maintain the existing support for WKT-encoded geometries in netCDF if there is now a more standard replacement ? The WKT support was perhaps a bit more generic than the new CF stuff that only supports "classic" simple geometries, whereas with WKT we could support exotic things like CircularString, CompoundCurve, CurvePolygon, etc.

[tomkralidis]

@rouault I think we're safe to move to the new convention (no backwards compat issues here).

[piyushrpt]

We are currently using OGR features to tag raster products with bounding box polygons. But we will be able to adopt to once support for new conventions is in place.

[wchen329]

It's good to see that no problems should be posed by the upgrade and that the upgrade will have further utility with developing new features. I've been currently getting showered with exams and other work at university, I hope next week I'll have something of importance to PR (as Spring Break finally arrives). Apologies for the delay.

Also, I understand that netCDFDataset::create( … ) [or something called within] is probably the entry point I want to use to make the important modifications, but I'm not entirely sure what is the most fitting manner to go implement the scanning of these extra attributes. It will be helpful to know about if there is some "catch all" routine where many of these "readers" such as (2306) FetchStandardParallels, some standard practices used for organizing these routines, and perhaps concurrency issues that may arise that aren't completely evident.

[dblodgett-usgs]

@wchen329 -- I've uploaded sample data (linked below)
Using that data, round trip the existing format with:

ogr2ogr -f netCDF old_ogr2ogr_states.nc states.shp
ogr2ogr -f "GeoJSON" states.geojson states.nc 

The file in here that is encoded in the current (CF-1.8) format was created with sf and ncdfgeom

states <- sf::read_sf("states.shp")
states <- ncdfgeom::write_geometry("ncdfgeom_states.nc", states)

states.tar.gz

[wchen329]

Hey thanks Dave (@dblodgett-usgs) ! Due to our conversation on Friday, I'll make sure to re-evaluate my approach some more before I tackle this again, but I think it will be easier using the approach you mentioned anyways.

[wchen329]

I've since written some utilities that take advantage of CF-1.8 encoding in NetCDF files and basically act as a high-level reader of CF-1.8 files with simple geometries. This is what I have wrote essentially (and has been tested and normal path of exec. valgrinded, though I have not pushed any testing utilities explicitly to GDAL)
https://github.com/wchen329/gdal/blob/master/gdal/frmts/netcdf/netcdfsg.h
https://github.com/wchen329/gdal/blob/master/gdal/frmts/netcdf/netcdfsg.cpp

At this point, I could kind of instead of messing with the existing grid-based and WKT-compliant GDAL NetCDF Driver, integrate this functionality into a new OGR driver. If one route is more or less useful for the GDAL project please let me know.

[rouault]

I'm not super keen in having a new netCDF driver. That would be confusing for users to have several netCDF drivers and would add unneeded technical debt to GDAL.
Based on @tomkralidis feedback, we could rework the existing vector functionality based on WKT (it would probably be great to probe the larger community by sending an email to the gdal-dev mailing list, although I would be surprised if anyone actually used that capability) and use CF-1.8 encoding instead.
For points, I believe that the support for featureType=point or featureType=profile could be left as it is currently.

[dblodgett-usgs]

Part of what @wchen329 is asking is if he should try to work with the existing code in the gdal-drivers folder that includes both gdal-raster functionality and WKT-based NetCDF functionality or add new code in an OGR driver folder that is specifically for NetCDF geometries?

The issue of deprecating or supporting the old WKT-based geometries approach is somewhat separate--or should be in my mind.

[rouault]

is if he should try to work with the existing code in the gdal-drivers folder that includes both gdal-raster functionality and WKT-based NetCDF functionality or add new code in an OGR driver folder that is specifically for NetCDF geometries?

Yes, I thought I had replied to that question in my above comment :-) My position would be to update the existing vector support in frmts/netcdf (the netcdflayer.cpp file specifically) to use CF-1.8 encoding in places where WKT is currently read/written. The WKT approach was my initial attempt at supporting non-point geometries since there was no other approach at the time where I implemented that. The existing driver has support for handling attribute fields, which hopefully can be kept as it is, supporting multiple layers in different files or in different netCDF groups, etc..

[dblodgett-usgs]

Ok. Thank you for the clarification.

That would be confusing for users

was confusing to me! Seems it would be more confusing to developers... no matter. The length of that file is a bit daunting, but we'll follow your guidance -- some more specific suggestions to make sure what @wchen329 implements would be a welcome PR would be helpful if you have some time to provide guidance.

[wchen329]

@dblodgett-usgs @rouault Thanks for the discussion and suggestions. From what has been said here it seems like introducing this sort of functionality to the existing driver perhaps as an additional condition check on the Conventions global attribute and allowing fall back on the WKT implementation for legacy datasets would be most beneficial. As @dblodgett-usgs mentioned some more specific suggestions on implementing this functionality could prove useful.

[dblodgett-usgs]

Sounds like a solid plan to me, @wchen329.

[dblodgett-usgs]

Tracking some test cases here. Sample cases and script to execute is included.
samples.zip

Some correspondence between @wchen329 and I for the record:

@dblodgett-usgs:

Looking mostly good.

1. For some reason, the spatial reference of the Yahara River shapefile is getting screwed up... maybe an X/Y mismatch? Open up the round tripped gpkg in QGIS and you'll see what I mean.
2. The states have a bunch of failures because of geometry type mismatches. Is this because multi polygon / polygon don't get handled right? You'll see the errors in the console and when you open the round trip gpkg in QGIS you'll see the missing states.

@wchen329:

The SRS feature is completely stolen from the CF-1.5 version of the driver so it'd be interesting to see where the issue lies. If it's a simple axis issue, it might be an issue with GIS.AxisOrderTraditional or whatever I used to correct the latitude/longitude issue. Or the CRS variable name might not be written correctly.
I'm not sure what the states issue could be. There is the thing where "Polygons" suddenly become "MultiPolygons" because of part node count and interior ring detection. Or maybe there's a more subtler thing occurring here. I'll make sure to look at these in a bit.

@dblodgett-usgs:

Regarding the Z/M thing. Could you just warn and drop the z/m content? Or add a flag to force drop the z/m? this data in particular doesn't actually have any Z or M coordinates.

@wchen329:

Sure we can do something that just shaves the M off of anything. 3D projected (Z) features types, however, are supported (so long as they don't have the M measure). I'll look into something in reducing the rank of the point to the "least non-lossy rank".

[wchen329]

@dblodgett-usgs thanks for posting this, I'll add what you have to the OP. Also just to clarify:

it might be an issue with GIS.AxisOrderTraditional

I mean to say here "an issue with how I use it" I don't suspect any issue within it itself.

[wchen329]

Okay, so I checked the states issue out, this time by converting to ESRI Shapefile. It looks like perhaps the enforcement of having the same Feature Type in the same layer is a little bit strict. Looking at it now. Perhaps the Yahara issue is related?

[dblodgett-usgs]

I don't follow what you mean by:

enforcement of having the same Feature Type in the same layer

[wchen329]

I don't follow what you mean by:

enforcement of having the same Feature Type in the same layer

So some drivers may choose to reduce some MultiPolygons to Polygons if they only had one polygon (which of course makes quite a bit of sense). Technically Polygons and MultiPolygons are considered separate wkbGeometry types. The restriction was put in place to prevent some potentially "less than ideal" -nlt usages in ogr2ogr, but The PR currently up removes the restriction (also because -nlt is not 100% guaranteed to make anything correct anyways, per the gdal docs). Anything mapping from non-POINT geometries into POINT will *most likely be quite lossy (as the original node counts is lost). But this is left as a responsibility to the user this time around.

* okay okay, maybe in some cases no

[rouault]

Closing this issue. CF 1.8 support is now there. If there are remaining issues, they should go into dedicated tickets