netcdf driver does not support irregular grids / GDAL GEOLOCATION arrays

[etourigny]

The CF standard allows for irregular 2D grids with "Two-dimensional coordinate variables" specified in the "coordinates" attribute (see CF sections 5.2 and 5.6).

CF projections currently unsupported by GDAL (#4285) can be supported via this mechanism, and be represented in GDAL via GEOLOCATION arrays. Current code does not support reading and writing GEOLOCATION arrays in netcdf files.

Attaching examples netcdf files in projected SRS with lat/lon values in "coordinates" attribute.

[etourigny]

In trunk (r23971) : support 2D GEOLOCATION arrays when a projected variable has coordinates attribute and supporting lon/at arrays .

This supports the coordinates "lon lat" attribute as described in CF section 5.6, only when longitude and latitude are present (as per cf-4.1 and 4.2) and the variable has a supporting grid_mapping (projected SRS) definition. If the variables are present, they are used as GEOLOCATION arrays (see example below).

The WRITE_LONLAT handling has changed slightly: default is NO, except if has geolocation ; with IF_NEEDED : write if has geoloc or is not CF projection. The latter is to ensure that if GDAL does not support the projection by mapping the CF attributes to OGC WKT, then the GEOLOCATION arrays are written to the CF file (and possibly read from input dataset).

If GEOLOCATION values are present in an input CF file, these will be used in export, instead of calculating with transform as previously.

Example (attached):

$ ncdump -h orog_RCM3.nc
netcdf orog_RCM3 {
dimensions:
	yc = 104 ;
	xc = 134 ;
variables:
	char Transverse_Mercator ;
		Transverse_Mercator:grid_mapping_name = "transverse_mercator" ;
		Transverse_Mercator:longitude_of_central_meridian = -97. ;
		Transverse_Mercator:latitude_of_projection_origin = 47.5 ;
		Transverse_Mercator:scale_factor_at_central_meridian = 1. ;
		Transverse_Mercator:false_easting = 3925000. ;
		Transverse_Mercator:false_northing = 3175000. ;
	double yc(yc) ;
		yc:long_name = "y-coordinate in Cartesian system" ;
		yc:standard_name = "projection_y_coordinate" ;
		yc:units = "m" ;
		yc:axis = "Y" ;
	double xc(xc) ;
		xc:long_name = "x-coordinate in Cartesian system" ;
		xc:standard_name = "projection_x_coordinate" ;
		xc:units = "m" ;
		xc:axis = "X" ;
	double lon(yc, xc) ;
		lon:units = "degrees_east" ;
		lon:long_name = "longitude" ;
		lon:standard_name = "longitude" ;
		lon:axis = "X" ;
	double lat(yc, xc) ;
		lat:units = "degrees_north" ;
		lat:long_name = "latitude" ;
		lat:standard_name = "latitude" ;
		lat:axis = "Y" ;
	float orog(yc, xc) ;
		orog:grid_mapping = "Transverse_Mercator" ;
		orog:standard_name = "surface_altitude" ;
		orog:long_name = "Surface Altitude" ;
		orog:units = "m" ;
		orog:coordinates = "lon lat" ;

// global attributes:
		:Conventions = "CF-1.0" ;
}

$ gdalinfo NETCDF:"orog_RCM3.nc":orog
Driver: netCDF/Network Common Data Format
Files: orog_RCM3.nc
Size is 134, 104
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",-97],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",3925000],
    PARAMETER["false_northing",3175000]]
Origin = (575000.000000000000000,5775000.000000000000000)
Pixel Size = (50000.000000000000000,-50000.000000000000000)
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
[...]
  yc#units=m
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]
  X_BAND=1
  X_DATASET=NETCDF:"orog_RCM3.nc":lon
  Y_BAND=1
  Y_DATASET=NETCDF:"orog_RCM3.nc":lat
Corner Coordinates:
Upper Left  (  575000.000, 5775000.000) ( 33d52' 4.42"W, 73d 9'59.02"N)
Lower Left  (  575000.000,  575000.000) (128d 1'21.79"W, 20d55'33.01"N)
Upper Right ( 7275000.000, 5775000.000) (160d 7'55.58"W, 73d 9'59.02"N)
Lower Right ( 7275000.000,  575000.000) ( 65d58'38.21"W, 20d55'33.01"N)
Center      ( 3925000.000, 3175000.000) ( 97d 0' 0.00"W, 47d30' 0.00"N)
Band 1 Block=134x1 Type=Float32, ColorInterp=Undefined
[...]

$ gdalwarp -geoloc -of netcdf NETCDF:"orog_RCM3.nc":orog  tmp2.nc
Creating output file that is 501P x 121L.
Processing input file NETCDF:orog_RCM3.nc:orog.
Using internal nodata values (eg. 1e+20) for image NETCDF:orog_RCM3.nc:orog.
0...10...20...30...40...50...60...70...80...90...100 - done.

Resulting lon/lat file is "similar" to original (TM re-projected to WGS84 lon/lat) although obviously not identical.

Due to bug in gdalwarp the extents are incorrect (resulting raster goes too far east), but the non-missing data pixels are correct.

Further testing with other datasets is needed.