Question regarding WCS DescribeCoverage ResponseCRSs

[ashleysommer]

I'm having trouble with QGIS interacting with our Thredds WCS endpoint. I'm not sure if its a bug in Thredds or QGIS or neither, or if I'm just interpreting it wrong.

Please see the NetCDF4 file served by Thredds 4.6.14 here:
http://esoil.io/thredds/wcs/SMIPSall/SMIPSv0.5.nc?service=WCS&version=1.0.0&request=GetCapabilities

Specifically the DescribeCoverage endpoint here: http://esoil.io/thredds/wcs/SMIPSall/SMIPSv0.5.nc?service=WCS&version=1.0.0&request=DescribeCoverage

The ResponseCRSs field contains:
EPSG:0 [Latitude_Longitude]

When QGIS queries the supportedCRSs from the DescribeCoverage endpoint, it sees the "EPSG:0" entry and interprets as an invalid CRS, so it ignores it.

When requesting data via a GetCoverage request, QGIS instead uses &responsecrs=OGC:CRS84 which of course causes Thredds to throw this error:
Response CRS [OGC:CRS84] not the supported CRS [EPSG:0 [Latitude_Longitude]]
If I substitute that with "EPSG:0 [Latitude_Longitude]" I get a valid response from WCS.

I feel like I'm missing something here. Where does the EPSG:0 [Latitude_Longitude] entry come from? Is it something I can change in my netcdf4 file? I have included a valid crs variable in the netcdf4 file, with correct spatialref metadata and CRS wkt entry.

Should QGIS interpret it as a valid CRS? I've opened a issue on the QGIS issue tracker about this here: qgis/QGIS#34569

This issue is related: qgis/QGIS#28864
But it seems to actually looks to be the opposite, it claims that when EPSG: 0 is found in the responseCRSs it should be ignored and fallback to the default QGIS value of OGC:CRS84.

[ashleysommer]

Ok, looks like its not something I'm doing wrong with our NetCDF4 file, even the Thredds example/sample WCS endpoint from the documentation has the same issue: https://thredds.ucar.edu/thredds/wcs/galeon/testdata/striped.nc?request=DescribeCoverage&version=1.0.0&service=WCS&coverage=ta

[ashleysommer]

I've done some digging, looks like this weird coordinate system does come from the netcdf4 file.

I am looking through the Thredds WCS code and found where EPSG: 0 is defined here:

thredds/cdm/src/main/java/ucar/unidata/geoloc/ogc/EPSG_OGC_CF_Helper.java

Line 66 in e96a745

Latitude_Longitude(0, "", "latitude_longitude"),

It mentions "Added in CF-1.2".
It is common (by convention) to add a variable named "crs" in the netcdf file, and add attributes describing that crs. One thing you can do is add "crs.grid_mapping_name=latitude_longitude" to specify that the grid used internally to the netcdf file is actually in latitude/longitude.

See example 5.8 in the CF-1.2 spec:
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.2/build/cf-conventions.html#coordinate-system

So when a netcdf4 file has "crs.grid_mapping_name=latitude_longitude" it appears that Thredds specifies the WCS responsecrs "EPSG:0 [Latitude_Longitude]"

FWIW, all of the relevant metadata in my netcdf4 file is the following:

Conventions=CF-1.6
spatial_ref='GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
crs.grid_mapping_name='latitude_longitude'
crs.long_name='Latitude Longitude'
crs.longitude_of_prime_meridian=0.0
crs.semi_major_axis=6378137.0
crs.inverse_flattening=298.257223563
crs.spatial_ref = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
crs.crs_wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'

[lesserwhirls]

Thank you for opening an issue @ashleysommer! I'm not super familiar with our WCS service (tagging @ethanrd), and the code is pretty old at this point, but I am almost certain we are handling this incorrectly. I could be wrong here, but in order for the TDS WCS to return the proper <responseCRS> code, we'd need to do one of two things:

1. Actually map the crs as defined by the CF parameters into a proper EPSG value
2. Look for the WKT definition in the CRS variable (crs.crs_wkt), trust that the data provider knows what they are doing, and map that to the proper EPSG code (probably more correct than CF → EPSG).
3. Provide a mechanism in the configuration catalog to allow the data provider the ability to associate the dataset with a specific EPSG code, and use that if it is set.

[ashleysommer]

Hi @lesserwhirls Thanks for your response.

I think there might be some of that logic already happening in the WCS code.
See, the <CoverageOffering>/<lonLatEnvelope.srs_name> has the correct SRS code to match the CRS given in the CF metadata.
Also the <CoverageOffering>/<domainSet>/<spatialDomain>/<EnvelopeWithTimePeriod> has correct code in SRS name.

I see the contents of <spatialDomain>/<RectifiedGrid> seem correct (if I understand correctly what that represents), however we see our friend EPSG: 0 on RectifiedGrid.srsName which I'm not sure is correct. The SRS on the RectifiedGrid should be some kind of x/y lookup.

Finally, the <supportedCRSs>/<requestCRSs> seem to have the correct value as per the netcdf CF metadata too.

So I think there's just a logic breakdown somewhere around determining the SRS for the <spatialDomain>/<RectifiedGrid> and somehow that makes its way into the <responseCRSs> entry.

I admit most of this is way over my head, so I might be completely wrong on all of that. And I don't really understand the difference between requestCRSs and responseCRSs in this case, or why they would be different.