NCSS return confusing for lat/lon data

[dopplershift]

For a straightforward request for data from models on lat/lon grids (say: http:https://thredds.ucar.edu/thredds/ncss/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_20180302_0000.grib2?var=Temperature_isobaric&disableLLSubset=on&disableProjSubset=on&horizStride=1&time_start=2018-03-02T00%3A00%3A00Z&time_end=2018-03-18T00%3A00%3A00Z&timeStride=1&vertCoord=&accept=netcdf4), the return doesn't seem to be CF-compliant:

netcdf GFS_Global_0p5deg_20180302_0000.grib2 {
dimensions:
	time1 = 93 ;
	isobaric = 31 ;
	lat = 361 ;
	lon = 720 ;
variables:
	float Temperature_isobaric(time1, isobaric, lat, lon) ;
		Temperature_isobaric:long_name = "Temperature @ Isobaric surface" ;
		Temperature_isobaric:units = "K" ;
		Temperature_isobaric:abbreviation = "TMP" ;
		Temperature_isobaric:missing_value = NaNf ;
		Temperature_isobaric:grid_mapping = "LatLon_Projection" ;
		Temperature_isobaric:coordinates = "time1 isobaric lat lon " ;
		Temperature_isobaric:Grib_Variable_Id = "VAR_0-0-0_L100" ;
		Temperature_isobaric:Grib2_Parameter = 0, 0, 0 ;
		Temperature_isobaric:Grib2_Parameter_Discipline = "Meteorological products" ;
		Temperature_isobaric:Grib2_Parameter_Category = "Temperature" ;
		Temperature_isobaric:Grib2_Parameter_Name = "Temperature" ;
		Temperature_isobaric:Grib2_Level_Type = 100 ;
		Temperature_isobaric:Grib2_Level_Desc = "Isobaric surface" ;
		Temperature_isobaric:Grib2_Generating_Process_Type = "Forecast" ;
	double time1(time1) ;
		time1:units = "Hour since 2018-03-02T00:00:00Z" ;
		time1:standard_name = "time" ;
		time1:long_name = "GRIB forecast or observation time" ;
		time1:calendar = "proleptic_gregorian" ;
		time1:_CoordinateAxisType = "Time" ;
	float isobaric(isobaric) ;
		isobaric:units = "Pa" ;
		isobaric:long_name = "Isobaric surface" ;
		isobaric:positive = "down" ;
		isobaric:Grib_level_type = 100 ;
		isobaric:_CoordinateAxisType = "Pressure" ;
		isobaric:_CoordinateZisPositive = "down" ;
	float lat(lat) ;
		lat:units = "degrees_north" ;
		lat:_CoordinateAxisType = "Lat" ;
		lat:standard_name = "latitude" ;
	float lon(lon) ;
		lon:units = "degrees_east" ;
		lon:_CoordinateAxisType = "Lon" ;
		lon:standard_name = "longitude" ;

// global attributes:
		:Originating_or_generating_Center = "US National Weather Service, National Centres for Environmental Prediction (NCEP)" ;
		:Originating_or_generating_Subcenter = "0" ;
		:GRIB_table_version = "2,1" ;
		:Type_of_generating_process = "Forecast" ;
		:Analysis_or_forecast_generating_process_identifier_defined_by_originating_centre = "Analysis from GFS (Global Forecast System)" ;
		:Conventions = "CF-1.6" ;
		:history = "Read using CDM IOSP GribCollection v3" ;
		:featureType = "GRID" ;
		:History = "Translated to CF-1.0 Conventions by Netcdf-Java CDM (CFGridWriter2)\nOriginal Dataset = /data/ldm/pub/native/grid/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_20180302_0000.grib2.ncx3#LatLon_361X720-p25S-180p0E; Translation Date = 2018-03-02T07:05:42.124Z" ;
		:geospatial_lat_min = -90. ;
		:geospatial_lat_max = 90. ;
		:geospatial_lon_min = 0. ;
		:geospatial_lon_max = 359.5 ;
}

Specifically, the grid_mapping attribute references a LatLon_Projection variable, which is not included in the file. While it's not necessary, it would probably be good to include this as a latitude_longitude projection so that the assumed radius of the earth can be included.

I have this this on 0.5 GFS and WW3 output. A netCDF file generated a year ago (for the 0.25 GFS) doesn't seem to have this problem--but only because there isn't a grid_mapping attribute. I'm not sure that's actually better.

[cwardgar]

This issue is related to #450 and #337: NCSS NetCDF output files only include selected data variables and the variables that comprise their coordinate systems. NCSS is not smart enough to examine attributes to find other variables that should be included. Auxiliary coordinate variables (i.e. variables listed in the coordinates attribute) are particularly tricky, because they would need to be subset in the same way that the primary coordinate variables are.

[dopplershift]

While I don't doubt the sentiment of your statement, the fact that this problem doesn't occur with e.g. Lambert Conformal output from HRRR on TDS would seems to indicate that the problem is specific to lat/lon coordinate systems, and not a more general NCSS issue.

[cwardgar]

Hmm, you're right. First, in addition to the variables I mentioned above, NCSS also includes the coordinate transform variables. However, when a GribCoverageDataset creates its CoverageCollections, transforms of type LatLon_Projection are explicitly excluded. So, as is often the case, the issue is upstream of NCSS.

Does anyone know why this decision was made? Perhaps John forgot that there are some cases where including a LatLon_Projection transform is actually necessary, e.g. when using a non-standard Earth radius? The fix is simple, but I'd like to understand the reasoning a bit better before I make the change and potentially break a bunch of stuff.

[dopplershift]

I'm guessing that's right. CF doesn't require it, but it also doesn't mind if you do include it. My preference is to include it always, because it simplifies my code (though I'll probably have to deal with lat/lon data without it from other source); if we leave it out, then we just need to make sure that it doesn't list a grid_mapping.

[lesserwhirls]

No idea why LatLon_Projection is excluded. Given that CoverageCollecion is so new, changing this behavior won't break anyone's code, and including it seems critical to actually using the data in a meaningful way. For example, without it, how is an end user going to compute spatial derivatives?

According to CF:

This grid mapping [grid_mapping_name = latitude_longitude] defines the canonical 2D geographical coordinate system based upon latitude and longitude coordinates on a spherical Earth. It is included so that the figure of the Earth can be described.

I'd say we should always include it for maximum clarity.