WMS with native grid

[PBrockmann]

Is it possible that a Thredds server could deliver a WMS from a netCDF (CF compliant) that has only only projected coordinates (here they are in a North Polar Stereographic projection - EPSG:32661) ?

I have a test file available from : https://thredds-su.ipsl.fr/thredds/catalog/ipsl_thredds/brocksce/tmp/catalog.html?dataset=DatasetScanIPSLTHREDDS/brocksce/tmp/g10.nc

But I get for the WMS : "No grids found in underlying NetCDF dataset"

The netCDF header looks like this:

netcdf GRISLI {
dimensions:
        time = UNLIMITED ; // (10 currently)
        y = 241 ;
        x = 241 ;
variables:
        double H(time, y, x) ;
                H:long_name = "ice_thickness" ;
                H:units = "m" ;
                H:coordinates = "x y" ;
        double time(time) ;
                time:long_name = "model_time" ;
                time:standard_name = "time" ;
                time:units = "days since 1900-1-1" ;
                time:calendar = "noleap" ;
        double x(x) ;
                x:long_name = "cartesian_x_coordinate" ;
                x:standard_name = "projection_x_coordinate" ;
                x:units = "m" ;
                x:axis = "x" ;
        double y(y) ;
                y:long_name = "cartesian_y_coordinate" ;
                y:standard_name = "projection_y_coordinate" ;
                y:units = "m" ;
                y:axis = "y" ;

// global attributes:
                :Conventions = "CF-1.8" ;

Here is a plot of the variable H from a cartopy python script:

[guygriffiths]

I believe that the issue here is that your data is not CF-compliant. From the CF standards: "When the horizontal spatial coordinate variables are not longitude and latitude, it is required that further information is provided to geo-locate the horizontal position. A grid mapping variable provides this information.". Currently there's nothing in your NetCDF file header that informs us what projection you're using.

The usual thing to do is have a pair of 2d lat/lon co-ordinates, but I assume that's not possible? I know there are other ways of defining the projection with an additional variable as well, but I don't have much experience with this.

What I can tell you with 100% certainty is that the message "No grids found in underlying NetCDF dataset" is coming directly from the Unidata NetCDF libraries (https://www.unidata.ucar.edu/software/netcdf-java/) which read the data. So ultimately there's nothing I can change in the code which will fix it (although I'm happy to help try to fix your data to solve the issue).

[ghansham]

Dear Sir

As it looks from ncml, the data in the file is map projected. But the projection information is missing in this file. See grid_mapping in the CF conventions document. For different map projections supported in CF conventions see this link:

https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/apf.html

Also see section 5.6 in the following link for how map projection information is packed in CF complaint netcdf files:

https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/cf-conventions.html#coordinate-system

Regards
Ghansham

[PBrockmann]

Could you point me real cases like the one from "Example 5.12. British National Grid + Newlyn Datum in CRS WKT format" ?
This is exactly what I would like to serve from our thredds server.

Indeed my netcdf test file has missing metadata "grip_mapping" attribute and crs variable (what value for this int variable BTW). I really miss real cases examples.

Will the WMS be available in this case, in other words will ncWMS will be able to digest and process those metadata ?

I have set a test file with the grid_mapping attribute and the a crs variable:
https://thredds-su.ipsl.fr/thredds/dodsC/ipsl_thredds/brocksce/tmp/g11.nc.html

ncWMS will be able to process those metadata and deliver a WMS ?

[ghansham]

Dear Sir

You are on the right path.

Actually I tried to repair the g10.nc using ncml approach:

https://docs.unidata.ucar.edu/thredds/ncml/2.2/ncml_cookbook.html

I have picked up the projection parameters from:

https://spatialreference.org/ref/epsg/32661/html/

In case the values of those projection parameters that you are using are different, you can modify the g10_repair.ncml.

And CF complaint projection information for Polar Stereographic has been picked up from:

https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/apf.html

See parameters in polar stereographic section. I have supplied them in the g10_repair.ncml

g10_repair.ncml.txt

Please note that, the false_easting and false_northing are in Kms which is not mentioned in CF convention document.
And "semi_major_axis" and "semi_minor_axis" have been added to specify the ellipsoidal shape of the earth. If they are not mentioned, then earth is considered spherical.

Please rename g10_repair.ncml.txt to g10_repair.ncml. Actually github does not support .ncml files to be uploaded.

THREDDS supports both netcdf and ncml files.

And also the ncml that I provided can be directly injected into thredds for repairing the datasets.
https://docs.unidata.ucar.edu/tds/current/userguide/using_ncml_in_the_tds.html

Also see the snapshot of the image in Unidata IDV (which uses NetCDF Java library, the same that is used by thredds and ncwms for reading netcdf files).

There is a very good desktop tool available at unidata (toolsUI):
https://downloads.unidata.ucar.edu/netcdf-java

There is good documentation for using toolsUI available here:

https://docs.unidata.ucar.edu/netcdf-java/current/userguide/toolsui_ref.html

You can go to the Feature Types Tab, within that go to Grids Tab.
You can check for the grids available by opening the desired file.

Let me know if you are facing any issue.

regards.
Ghansham

[PBrockmann]

Thanks very much for helping Ghansham. Your return experience is very valuable.

I have corrected directly my test file from a CDL file editing but I note the use of the ncml additionnal file.

Here are the new header:

netcdf g12 {
dimensions:
        time = UNLIMITED ; // (10 currently)
        y = 241 ;
        x = 241 ;
variables:
        double H(time, y, x) ;
                H:long_name = "ice_thickness (cf,Searise)" ;
                H:units = "m" ;
                H:coordinates = "x y" ;
                H:grid_mapping = "crs" ;
        double time(time) ;
                time:long_name = "model_time" ;
                time:standard_name = "time" ;
                time:units = "days since 1900-1-1" ;
                time:calendar = "noleap" ;
        double x(x) ;
                x:long_name = "cartesian_x_coordinate" ;
                x:standard_name = "projection_x_coordinate" ;
                x:units = "m" ;
                x:axis = "x" ;
        double y(y) ;
                y:long_name = "cartesian_y_coordinate" ;
                y:standard_name = "projection_y_coordinate" ;
                y:units = "m" ;
                y:axis = "y" ;
        int crs ;
                crs:grid_mapping_name = "Polar Stereographic North" ;
                crs:proj4_params = "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs +type=crs" ;
                crs:EPSG_code = "EPSG:32661" ;
                crs:straight_vertical_longitude_from_pole = 0 ;
                crs:latitude_of_projection_origin = 90 ;
                crs:scale_factor_at_projection_origin = 0.994 ;
                crs:false_easting = 2000 ;
                crs:false_northing = 2000 ;
                crs:semi_major_axis = 6378137 ;
                crs:semi_minor_axis = 6356752.5 ;

// global attributes:
                :Conventions = "CF-1.8" ;

But the WMS is still not available: "No grids found in underlying NetCDF dataset"

https://thredds-su.ipsl.fr/thredds/wms/ipsl_thredds/brocksce/tmp/g12.nc?service=WMS&version=1.3.0&request=GetCapabilities

What do I have missed ?

[guygriffiths]

You've missed the fact that you can't just pick an arbitrary name for the grid mapping, you need to choose one from this list: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/apf.html

Specifically, you need to choose polar_stereographic, as Ghansham suggested.

However, in trying this out, I found a bug in ncWMS which means this doesn't display correctly (it incorrectly calculates the bounding box of the data, causing it to not display above latitudes of 39). I'll look into this at some point next week, but it's unlikely to work properly using THREDDS. Give it a try though and let me know how you get on

[ghansham]

Dear Sir

@guygriffiths sir is absolutely correct. The "grid_mapping_name" attribute of the "crs" variable (in this case) must be "polar_stereographic". Rest everything seems correct. Please refer to the ncml that I shared in my previous post.

Regards
Ghansham

[PBrockmann]

Thanks again on your interest on this topic.
I try to understand what are the metadata that are really compulsory.
I have set the minimal attributes to get the WMS available.

Here they are:

netcdf g14 {
dimensions:
        time = UNLIMITED ; // (10 currently)
        y = 241 ;
        x = 241 ;
variables:
        double H(time, y, x) ;
                H:long_name = "ice_thickness (cf,Searise)" ;
                H:units = "m" ;
                H:coordinates = "x y" ;
                H:grid_mapping = "crs" ;
        double time(time) ;
                time:long_name = "model_time" ;
                time:standard_name = "time" ;
                time:units = "days since 1900-1-1" ;
                time:calendar = "noleap" ;
        double x(x) ;
                x:long_name = "cartesian_x_coordinate" ;
                x:standard_name = "projection_x_coordinate" ;
                x:units = "m" ;
                x:axis = "x" ;
        double y(y) ;
                y:long_name = "cartesian_y_coordinate" ;
                y:standard_name = "projection_y_coordinate" ;
                y:units = "m" ;
                y:axis = "y" ;
        int crs ;
                crs:grid_mapping_name = "polar_stereographic" ;
                crs:EPSG_code = "EPSG:32661" ;
                crs:straight_vertical_longitude_from_pole = 0 ;
                crs:latitude_of_projection_origin = 90 ;
                crs:scale_factor_at_projection_origin = 0.994 ;

// global attributes:
                :Conventions = "CF-1.8" ;

The WMS is available from https://thredds-su.ipsl.fr/thredds/wms/ipsl_thredds/brocksce/tmp/g15.nc?service=WMS&version=1.3.0&request=GetCapabilities

The representation I get with a simple leaflet code is strange:

See:

<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<style>
  .map {
    height: 600px;
    width: 800px;
  }
</style>

<script src="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.9.4/leaflet.js"></script>
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.9.4/leaflet.css" type="text/css">

<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.9.2/proj4.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4leaflet/1.0.2/proj4leaflet.js"></script>

</head>

<body>
<h2>My Map</h2>
<div id="map" class="map"></div>
<script type="text/javascript">

//============================================
var crs1 = new L.Proj.CRS('EPSG:32661',
        '+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs +type=crs',
        {
                resolutions: [32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8],
        }
);

map = new L.Map('map', {
        crs: crs1,
        center: [90, 0],                        // lat, lon
        zoom: 0,
});

//============================================

var ressource = "https://thredds-su.ipsl.fr/thredds/wms/ipsl_thredds/brocksce/tmp/g14.nc";
var wmsLayer = L.tileLayer.wms(
    ressource + '?', {
        layers: 'H',
        elevation: 0,
        time: '1831-03-27T00:00:00.000Z',
        colorscalerange: '0,4150',
        numcolorbands: 40,
        styles: 'raster/psu-plasma',
        abovemaxcolor: 'extend',
        belowmincolor: 'extend',
        format: 'image/png',
        opacity: 100
     }
);
wmsLayer.addTo(map);

//============================================

</script>
</body>
</html>

[guygriffiths]

I found a bug in ncWMS which means this doesn't display correctly (it incorrectly calculates the bounding box of the data, causing it to not display above latitudes of 39)

The representation I get with a simple leaflet code is strange:

The representation you're getting with leaflet is cutting everything above the latitude of 39 off, because ncWMS thinks that's the limit of the dataset. The reason for that is that the Unidata library is incorrectly reporting back the lat/lon bounding box of the projected co-ordinate system.

My best guess is that they are just checking the edges of the projected grid to determine the maximum available latitude. If we look at your leaflet plot, you can see that the centre of right-hand edge is the highest latitude point on the image perimeter. But of course in a polar projection, the actual highest latitude point is somewhere in the middle of the image.

It's not a big deal to fix, and I should be able to make a new release of ncWMS next week. That won't update THREDDS as far as I'm aware, but you may be able to simply drop the new library in place of the old one on your Tomcat installation.

[ghansham]

Dear sir

You should check the repair ncml that I sent for minimal amount of metadata required to correct your file. I think you haven't checked that. In the latest metadata that you posted has skipped semi_major and semi_minor axes. And EPSG_code is additional (not required).

Regards
Ghansham

[PBrockmann]

Hi Ghansham,
I have tried to see what are the minimal metadate that are needed by ncWMS to deliver the service.
Where is it documented ?
Compare g13.nc and g14.nc. The use of g13.nc with my leaflet test script and the problem is still there.
https://thredds-su.ipsl.fr/thredds/dodsC/ipsl_thredds/brocksce/tmp/g13.nc.html

[guygriffiths]

@PBrockmann - The data which is not working in leaflet, due to the ncWMS bug, may only not be buggy because leaflet works exclusively in EPSG:3857 and uses lat-lon. I suspect that if you actually make an EPSG:32661 request to ncWMS directly it may work. The fact that ncWMS is getting the incorrect lat-lon bounding box shouldn't come into play if we never translate into lat-lon (I'm not 100% certain on this, but I would definitely recommend trying it).

[PBrockmann]

For info the bug also appears in EPSG:32761

See:
https://thredds-su.ipsl.fr/thredds/catalog/ipsl_thredds/brocksce/tmp/catalog.html?dataset=DatasetScanIPSLTHREDDS/brocksce/tmp/GRISLI_Antarctica_02.nc

//============================================
var crs1 = new L.Proj.CRS('EPSG:32761',
        '+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs +type=crs',
        {
                resolutions: [32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8],
        }
);

map = new L.Map('map', {
        crs: crs1,
        center: [-90, 0],                        // lat, lon
        zoom: 0,
});

//============================================

var ressource = "https://thredds-su.ipsl.fr/thredds/wms/ipsl_thredds/brocksce/tmp/GRISLI_Antarctica_02.nc";
var wmsLayer = L.tileLayer.wms(
    ressource + '?', {
        layers: 'H',
        elevation: 0,
        time: '1900-01-01T00:00:00.000Z',
        colorscalerange: '0,4150',
        numcolorbands: 40,
        styles: 'raster/psu-plasma',
        abovemaxcolor: 'extend',
        belowmincolor: 'extend',
        format: 'image/png',
        opacity: 100
     }
);
wmsLayer.addTo(map);

Indeed I could also describe the grid with latitudes, longitudes expressed in EPSG:4326 but my purpose is to explore if I can
display map with WMS at the native grid.

[guygriffiths]

No, I think not until I publish a new version of ncWMS with the bug fixed, which will be next week sometime.

[ghansham]

Sir

You have to follow CF conventions document. That what we follow.

Regards
Ghansham

[PBrockmann]

Any new code to test on this issue ?

[sbjartmar]

No, I think not until I publish a new version of ncWMS with the bug fixed, which will be next week sometime.

Did you release a version somewhere with this fix?

[guygriffiths]

No, apologies but I haven't had time to work much on this. There is a snapshot release of the edal-cdm library which fixes this issue though. You can download it from https://github.com/Reading-eScience-Centre/edal-java/releases/tag/edal-1.5.3-SNAPSHOT and drop it in in place of edal-1.5.2.jar in the WEB-INF directory.

[PBrockmann]

Many thanks.

I have tested and indeed the map appears but with gap. A "closing line" appears whereas there should not have any since there is no reprojection needed.

I have used OpenLayers and the g14.nc file available from
https://thredds-su.ipsl.fr/thredds/catalog/ipsl_thredds/brocksce/tmp/catalog.html?dataset=DatasetScanIPSLTHREDDS/brocksce/tmp/g14.nc

Code used:
https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/map_g14.html

That produces

Other point, is the patch possible to be applied for the ncWMS code distributed from a Thredds server (ie code coming from the thredds.war) ?

To test your correction I have copied edal-cdm-1.5.3-SNAPSHOT.jar to /usr/local/tomcat/webapps/ncWMS2/WEB-INF/lib/edal-cdm-1.5.2.jar.
What about the /usr/local/tomcat/webapps//thredds/WEB-INF/lib/edal-cdm-1.5.0.5-SNAPSHOT.jar ?

[PBrockmann]

Tested with OpenLayers.

Code used:
https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/map_g14_ol.html (install OpenLayers with npm).

Same closing problem.

[PBrockmann]

Any news on this issue ?

[PBrockmann]

Am I the only one to be interested to use WMS at native grids ?

[guygriffiths]

Hi, sorry but I'm not really actively working on this at the moment. I'm happy to review any PRs which fix the issue, but I don't have time to look into it myself.

[PBrockmann]

Ok could you tell me where to look at on this issue nearly corrected ? The last problem that remains is about non closing grid at the change date line that left some parts unplotted.
My purpose is to display a grid at its native definition. The ice land model works at the specified projection.

[guygriffiths]

Sure, I think the entry-point to start looking is the drawIntoImage() method defined in https://github.com/Reading-eScience-Centre/edal-java/blob/master/graphics/src/main/java/uk/ac/rdg/resc/edal/graphics/style/ImageLayer.java

I would stick a breakpoint there and step through from that point to get to the issue.