Unstructured grids - UGRID from a TDS

[PBrockmann]

I have read this proposal and I am wondering if the project https://github.com/ioos/thredds-ugrid has been finally merged into the current Thredds code.

This proposal was a really good initiative.
https://github.com/ioos/thredds-ugrid/blob/master/documents/proposal/2018-0196%20Proposal_Narrative_redacted.pdf
How did this project ended ?

I am looking for an infrastructure able to deliver a WMS for UGRID ressources.
Let me know.

[haileyajohnson]

I believe this work is more or less complete but still needs testing... @lesserwhirls can probably provide more information

[PBrockmann]

I would be really interested to test it.
There are a couple of files from https://geoport.whoi.edu/ but can provide output from our icosahedral hydrostatic model.

[PBrockmann]

Any news ? Having a calendar in mind would be kind.

[lesserwhirls]

The general work is complete, and FMRCs can be created using UGRID datasets. The work is not available in the current release (v5.4), but will be in the next release (and should be available in the current v5.5 SNAPSHOTS, available at https://downloads.unidata.ucar.edu/tds/). The WCS service will not work with UGRID datasets, and NCSS currently does not work with UGRID datasets (but will be with a bit more work). An example of a UGRID FMRC can be found here:

tds/tds/src/test/content/thredds/catalogFmrc.xml

Lines 78 to 87 in dda75f8

	<featureCollection featureType="FMRC" name="USF_WSF_nc" path="testUsfWsfUgridFmrc">
	<metadata inherited="true">
	<dataType>Grid</dataType>
	<dataFormat>NetCDF</dataFormat>
	<serviceName>ugrid</serviceName>
	</metadata>
	<collection spec="${cdmUnitTest}/conventions/ugrid/USF_WSF/usf_fvcom_#yyyy_MM_dd#.nc$" />
	<update startup="true" trigger="allow"/>
	<fmrcConfig regularize="true" datasetTypes="TwoD Best Files Runs ConstantForecasts ConstantOffsets"/>
	</featureCollection>

The definition of the "ugrid" service is:

tds/tds/src/test/content/thredds/catalog.xml

Lines 24 to 32 in dda75f8

	<service name="ugrid" base="" serviceType="compound">
	<service name="odap" serviceType="OpenDAP" base="/thredds/dodsC/"/>
	<service name="http" serviceType="HTTPServer" base="/thredds/fileServer/"/>
	<service name="wms" serviceType="WMS" base="/thredds/wms/"/>
	<service name="cdmremote" serviceType="CdmRemote" base="/thredds/cdmremote/"/>
	<service name="iso" serviceType="ISO" base="/thredds/iso/"/>
	<service name="ncml" serviceType="NCML" base="/thredds/ncml/"/>
	<service name="uddc" serviceType="UDDC" base="/thredds/uddc/"/>
	</service>

[PBrockmann]

Very good.
Is the WMS will be available for UGRID datasets ?

[lesserwhirls]

Very good.

Is the WMS will be available for UGRID datasets ?

Yes.

[PBrockmann]

Impatient to test it.
Do you have a TDS available for testing ?

[lesserwhirls]

I do not - I am no longer an employee at Unidata, and I do not run a TDS from home. Docker might be the easiest way to test.

[haileyajohnson]

Impatient to test it. Do you have a TDS available for testing ?

You can set up your own instance by downloading the latest snapshot here, or you can see a running instance of the latest snapshot at thredds-test.unidata.ucar.edu - I don't know if there are any UGRID datasets on it currently though.

[PBrockmann]

Ok thank you. I will investigate. Regards

[PBrockmann]

Ok tested. Indeed it works with the tested UGRID files I had.

Just how to debug when the GetCapabilities answers
Cause: No grids found in underlying NetCDF dataset

[PBrockmann]

It would be interesting to reopen this issue because I think there are some precisions to make.
Perhaps a GitHub discussion would be better.

I have tried a UGRID found from https://psyplot.github.io/examples/maps/example_ugrid.html
The ugrid_demo.nc. The WMS GetCapacities response seems ok but there are in fact no Layers available.

http://0.0.0.0/thredds/wms/testAll/ugrid_demo.nc?service=WMS&version=1.3.0&request=GetCapabilities

<Layer>
<Title>netCDF output from StormFlash2d</Title>
</Layer>

Is this file uncorrect from the UGRID convention ?

[lesserwhirls]

The code that handles interpreting the UGRID convention is located in the tds-ugrid subprojtect of the TDS codebase (see here). Also, I just noticed that you are dealing with hexahedron cells - this implementation only handles triangular meshes, so that's why no layers are found. I'm surprised there are no error messages in the logs about the face node dimensionality (I'll have to check that out).

[PBrockmann]

Ahhh. Only triangles. I understand better.

netcdf selfe_ugrid_select1 {
dimensions:
        face = 3 ;
        max_face_nodes = 3 ;
        node = 5 ;
variables:
        int Mesh ;
                Mesh:cf_role = "mesh_topology" ;
                Mesh:long_name = "Topology data of 2D unstructured mesh" ;
                Mesh:topology_dimension = 2 ;
                Mesh:node_coordinates = "node_lon node_lat" ;
                Mesh:face_node_connectivity = "face_nodes" ;
        int face_nodes(face, max_face_nodes) ;
                face_nodes:cf_role = "face_node_connectivity" ;
                face_nodes:long_name = "Maps every face to its corner nodes." ;
                face_nodes:start_index = 1 ;
        double node_lat(node) ;
                node_lat:standard_name = "latitude" ;
                node_lat:long_name = "Latitude of 2D mesh nodes" ;
                node_lat:units = "degrees_north" ;
        double node_lon(node) ;
                node_lon:standard_name = "longitude" ;
                node_lon:long_name = "Longitude of 2D mesh nodes" ;
                node_lon:units = "degrees_east" ;
        float surface_temp(node) ;
                surface_temp:coordinates = "node_lon node_lat" ;
                surface_temp:long_name = "Sea surface temperature" ;
                surface_temp:standard_name = "sea_surface_temperature" ;
                surface_temp:mesh = "Mesh" ;
                surface_temp:location = "node" ;
                surface_temp:units = "degree_Celsius" ;

// global attributes:
                :Conventions = "CF-1.6/UGRID-1.x" ;

data:

 Mesh = _ ;

 face_nodes =
  1, 4, 5,
  1, 5, 2,
  4, 1, 3;

 node_lat = 46.3698348999023, 46.3674011230469,
         46.3713722229004, 46.3680648803711, 46.3661003112793;

 node_lon = -124.104789733887, -124.100074768066,
         -124.110336303711, -124.110466003418, -124.105400085449;

 surface_temp = 200, 100, 0, 50, 80, 150;
}

Gives this.

But when I remove 2 nodes to have only a triangle, nothing appears.

netcdf selfe_ugrid_select1 {
dimensions:
        face = 1 ;
        max_face_nodes = 3 ;
        node = 3 ;
variables:
        int Mesh ;
                Mesh:cf_role = "mesh_topology" ;
                Mesh:long_name = "Topology data of 2D unstructured mesh" ;
                Mesh:topology_dimension = 2 ;
                Mesh:node_coordinates = "node_lon node_lat" ;
                Mesh:face_node_connectivity = "face_nodes" ;
        int face_nodes(face, max_face_nodes) ;
                face_nodes:cf_role = "face_node_connectivity" ;
                face_nodes:long_name = "Maps every face to its corner nodes." ;
                face_nodes:start_index = 1 ;
        double node_lat(node) ;
                node_lat:standard_name = "latitude" ;
                node_lat:long_name = "Latitude of 2D mesh nodes" ;
                node_lat:units = "degrees_north" ;
        double node_lon(node) ;
                node_lon:standard_name = "longitude" ;
                node_lon:long_name = "Longitude of 2D mesh nodes" ;
                node_lon:units = "degrees_east" ;
        float surface_temp(node) ;
                surface_temp:coordinates = "node_lon node_lat" ;
                surface_temp:long_name = "Sea surface temperature" ;
                surface_temp:standard_name = "sea_surface_temperature" ;
                surface_temp:mesh = "Mesh" ;
                surface_temp:location = "node" ;
                surface_temp:units = "degree_Celsius" ;

// global attributes:
                :Conventions = "CF-1.6/UGRID-1.x" ;

data:

 Mesh = _ ;

 face_nodes =
  1, 2, 3;

 node_lat = 46.3698348999023, 46.3674011230469,
         46.3713722229004;

 node_lon = -124.104789733887, -124.100074768066,
         -124.110336303711;

 surface_temp = 200, 100, 0;
}

Why ? And last question how to have value on faces, not on nodes ? An example would be welcomed.

[lesserwhirls]

In order for a variable to be recognized as a "gridded" type by netCDF-Java, the variable will need to have a time dimension and time coordinate variable. As for an example, check out the netCDF file and its associated NcML file here.

[PBrockmann]

Not sure to be agree. A variable can be time independant (area of cells, bathymetry, ...).
As a proof, the cdl selfe_ugrid_select1 I have put above is accepted and the WMS is available for mapping.

[lesserwhirls]

My apologies, you are correct - I was thinking of restrictions on a different TDS service.

[PBrockmann]

Here is what I get from this file that has no time dimension.

netcdf selfe_ugrid3 {
dimensions:
        face = 3 ;
        max_face_nodes = 3 ;
        node = 5 ;
variables:
        int Mesh ;
                Mesh:cf_role = "mesh_topology" ;
                Mesh:long_name = "Topology data of 2D unstructured mesh" ;
                Mesh:topology_dimension = 2 ;
                Mesh:node_coordinates = "node_lon node_lat" ;
                Mesh:face_node_connectivity = "face_nodes" ;
        int face_nodes(face, max_face_nodes) ;
                face_nodes:cf_role = "face_node_connectivity" ;
                face_nodes:long_name = "Maps every face to its corner nodes." ;
                face_nodes:start_index = 1 ;
        double node_lat(node) ;
                node_lat:standard_name = "latitude" ;
                node_lat:long_name = "Latitude of 2D mesh nodes" ;
                node_lat:units = "degrees_north" ;
        double node_lon(node) ;
                node_lon:standard_name = "longitude" ;
                node_lon:long_name = "Longitude of 2D mesh nodes" ;
                node_lon:units = "degrees_east" ;
        float surface_temp(node) ;
                surface_temp:coordinates = "node_lon node_lat" ;
                surface_temp:long_name = "Sea surface temperature" ;
                surface_temp:standard_name = "sea_surface_temperature" ;
                surface_temp:mesh = "Mesh" ;
                surface_temp:location = "node" ;
                surface_temp:units = "degree_Celsius" ;

// global attributes:
                :Conventions = "CF-1.6/UGRID-1.x" ;

data:

 Mesh = _ ;

 face_nodes =
  1, 4, 5,
  1, 2, 3,
  1, 3, 4;

 node_lat = 40, 50,
         50, 40, 30;

 node_lon = 0, 10,
         20, 30, 10;

 surface_temp = 200, 100, 0, 50, 150;

}

But I am not able to reduce to only 4 or 3 nodes. Do not understand why ?
I am also looking for testing face values (not on nodes).

[lesserwhirls]

I will need some free time on the weekend to dig into why exactly you are unable to reduce your example dataset (not sure when I will have the time though). In the meantime, have a look at this netCDF file, which works with the TDS and edal-java (provides WMS services) without needing any NcML modifications. This example file was produced by fvcom over the west Florida shelf and includes both node and face variables.

[PBrockmann]

Thank you Sean.
What I am missing is very simple examples.
Here working on how to represent 3 triangles with face values.
I am looking into your file. I should be ok with that. I will let you know.

[PBrockmann]

Nope. Even by inspecting the usf_fvcom_2019_01_28.nc that by the way do not respect the UGRID convention when I use the
ugrid-checker

$ ugrid-checker usf_fvcom_2019_01_28.nc
UGRID conformance checks complete.

List of checker messages :
  *** FAIL R202 : Mesh coordinate variable "lonc" within fvcom_mesh:face_coordinates has dimension "nele", but the parent mesh face dimension is "three".
  *** FAIL R202 : Mesh coordinate variable "latc" within fvcom_mesh:face_coordinates has dimension "nele", but the parent mesh face dimension is "three".
  *** FAIL R504 : Mesh data variable "hc" has location="cell", which is not one of "face", "edge" or "node".
  *** FAIL R509 : Mesh data variable "hc" has dimensions ('nele',), of which 0 are mesh dimensions, instead of 1.
  *** FAIL R509 : Mesh data variable "u" has dimensions ('time', 'siglay', 'nele'), of which 0 are mesh dimensions, instead of 1.
  *** FAIL R509 : Mesh data variable "v" has dimensions ('time', 'siglay', 'nele'), of which 0 are mesh dimensions, instead of 1.

Total of 6 problems logged :
  6 Rxxx requirement failures
  0 Axxx advisory recommendation warnings

Done.

So I don't get it. Here is a simple file with 3 triangles that I cannot plot with the WMS.
I get triangles but cut.

Here is the cdl.

netcdf ugrid9 {
dimensions:
        face = 3 ;
        node = 5 ;
        three = 3 ;
variables:
        int Mesh ;
                Mesh:cf_role = "mesh_topology" ;
                Mesh:long_name = "Topology data of 2D unstructured mesh" ;
                Mesh:topology_dimension = 2 ;
                Mesh:node_coordinates = "node_lat node_lon" ;
                Mesh:face_node_connectivity = "face_nodes" ;
                Mesh:face_coordinates = "lat_c lon_c" ;
        int face_nodes(three, face) ;
                face_nodes:cf_role = "face_node_connectivity" ;
                face_nodes:long_name = "Maps every face to its corner nodes." ;
                face_nodes:start_index = 1 ;
        double node_lat(node) ;
                node_lat:standard_name = "latitude" ;
                node_lat:long_name = "Latitude of 2D mesh nodes" ;
                node_lat:units = "degrees_north" ;
        double node_lon(node) ;
                node_lon:standard_name = "longitude" ;
                node_lon:long_name = "Longitude of 2D mesh nodes" ;
                node_lon:units = "degrees_east" ;
        double lat_c(face) ;
                lat_c:standard_name = "latitude" ;
                lat_c:long_name = "Latitude of 2D mesh nodes" ;
                lat_c:units = "degrees_north" ;
        double lon_c(face) ;
                lon_c:standard_name = "longitude" ;
                lon_c:long_name = "Longitude of 2D mesh nodes" ;
                lon_c:units = "degrees_east" ;
        float surface_temp(face) ;
                surface_temp:coordinates = "node_lon_c node_lat_c" ;
                surface_temp:long_name = "Sea surface temperature" ;
                surface_temp:standard_name = "sea_surface_temperature" ;
                surface_temp:mesh = "Mesh" ;
                surface_temp:location = "face" ;
                surface_temp:units = "degree_Celsius" ;

// global attributes:
                :Conventions = "CF-1.6/UGRID-1.0" ;

data:

 Mesh = _ ;

 face_nodes =
  1, 5, 4,
  1, 3, 2,
  1, 4, 3;

 node_lat = 40, 50, 50, 40, 30;
 node_lon = 0, 10, 20, 30, 10;
 lat_c = 35, 45, 45;
 lon_c = 20, 7.5, 15;

 surface_temp = 200, 100, 0;

}

I do not understand what is uncorrect. Do the centers of the faces need to be at the centroïds ?

[PBrockmann]

A simplified version of your file usf_fvcom_2019_01_28.nc gives https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/usf_fvcom_2019_01_28_extract2.nc (only one variable on faces).
and this map:

What is odd is that it also gives problems with the ugrid-checker (https://pypi.org/project/ugrid-checks/)

$ ugrid-checker usf_fvcom_2019_01_28_extract2.nc

UGRID conformance checks complete.

List of checker messages :
  *** FAIL R202 : Mesh coordinate variable "lonc" within fvcom_mesh:face_coordinates has dimension "nele", but the parent mesh face dimension is "three".
  *** FAIL R202 : Mesh coordinate variable "latc" within fvcom_mesh:face_coordinates has dimension "nele", but the parent mesh face dimension is "three".
  *** FAIL R509 : Mesh data variable "u" has dimensions ('nele',), of which 0 are mesh dimensions, instead of 1.

Total of 3 problems logged :
  3 Rxxx requirement failures
  0 Axxx advisory recommendation warnings

Done.

Error disappears if you change face_node_connectivity dimension order:
int nv(three, nele) to int nv(nele, three) but then the map is not correct.
The ugrid-checks tool may not be agreed to the reading from the ncWMS UGRID.
See also ugrid-conventions/ugrid-conventions#57

[PBrockmann]

What would really help me is to have two single triangle UGRID netcdf files. One with node values (3) and one a face value (1).

[PBrockmann]

@lesserwhirls you wrote: "this implementation only handles triangular meshes"
So hexaedron cannot be represented ?
That might be a problem for representing DYNAMICO output that uses an icosahedral grid.

[PBrockmann]

Any news on representing a single triangle with a face value ?
This simple representation would be help a lot to understand how TDS can represent a UGRID variable.

[lesserwhirls]

I've finally had some spare time to dig into this a bit. For a single triangle (three nodes, one face), the following cdl should generate an ok netCDF file:

netcdf bermuda_triangle {
dimensions:
        node = 3 ;
        face = 1 ;
        three = 3 ;
        time = 1 ;
variables:
        int fvcom_mesh ;
                fvcom_mesh:cf_role = "mesh_topology" ;
                fvcom_mesh:topology_dimension = 2 ;
                fvcom_mesh:node_coordinates = "lon lat" ;
                fvcom_mesh:face_node_connectivity = "nv" ;
        int time(time) ;
                time:units = "seconds since 1970-01-01" ;
                time:standard_name = "time" ;
        double lat(node) ;
                lat:units = "degrees_north" ;
        double latc(face) ;
                latc:units = "degrees_north" ;
        double lon(node) ;
                lon:units = "degrees_east" ;
        double lonc(face) ;
                lonc:units = "degrees_east" ;
        int nv(face, three) ;
                nv:long_name = "nodes surrounding element" ;
                nv:cf_role = "face_node_connectivity" ;
                nv:start_index = 0 ;
        double h(node) ;
                h:coordinates = "lat lon" ;
                h:mesh = "fvcom_mesh" ;
                h:location = "node" ;
                h:units = "m" ;
        float temp(time, node) ;
                temp:coordinates = "time lat lon" ;
                temp:mesh = "fvcom_mesh" ;
                temp:location = "node" ;
                temp:units = "degree_Celsius" ;
        float u(time, face) ;
                u:coordinates = "time latc lonc" ;
                u:mesh = "fvcom_mesh" ;
                u:location = "face" ;
                u:units = "m s-1" ;

// global attributes:
                :Conventions = "CF-1.X/UGRID-1.X" ;
data:

 fvcom_mesh = _ ;

 time = 1681674619 ;

 lat = 27.789, 18.466, 32.295 ;

 latc = 25.52 ;

 lon = -80.227, -66.106, -64.781 ;

 lonc = -70.37 ;

 nv =
  0, 1, 2 ;

 h = 9.1, 10.2, 11.3 ;

 temp =
  20.1, 21.2, 23.3 ;

 u =
  -0.0775 ;
}

However, if I am reading the edal-java code correctly, there is a bad assumption being made here:

https://github.com/Reading-eScience-Centre/edal-java/blob/e881f5deb8d0c219f305418ba07803e1546a1cb4/cdm/src/main/java/uk/ac/rdg/resc/edal/dataset/cdm/CdmGridDatasetFactory.java#L837-L859

For triangular meshes, the face_node_connecitivity variable should have dimensionality number_of_faces, 3. In the case of a single triangle (one face, each face with three edges), the edal-java code will incorrectly swap dimensions and throw errors. That's likely why you cannot go below four nodes (as two faces would trigger the incorrect behavior). I'd suggest filing a bug with the edal-java project.

[PBrockmann]

Thank you very much.
I will test this.

Is there any chance to have a hexahedron be represented ?

[lesserwhirls]

You would need to work with the edal-java project regarding hexahedron meshes in WMS, as their code would need to be extended.

[PBrockmann]

No sucess when I try to display the bermuda_triangle.
Nothing visible.

[lesserwhirls]

Right - there is a bug in edal-java that is preventing the single triangle case from working. See the last part of #345 (comment).

[PBrockmann]

Yes, it should be a good thing to describe this bug in an issue on the edal-java project
But there is already one about the subject - see Reading-eScience-Centre/edal-java#121, but still open since Oct 2019...