GDAL's NetCDF driver doesn't seem to read the spatial reference system as documented

[hiebert]

I'm trying to use GDAL (version 1.7.2 on Gentoo built with hdf5 and netcdf support) to read NetCDF files, but I'm a little unclear on how the NetCDF driver attempts to determine the georeferencing.

The driver page here: ​http://www.gdal.org/frmt_netcdf.html states that it will attempt to read the metadata "spatial_ref" for a Well Know Text definition of the spatial reference (I assume that it means that it will try to use a global attribute of the NetCDF file named spatial_ref).

Since this would be convenient for my application, I attempted to use GDAL this way without success:

• I added a global attribute named "spatial_ref" with the WKT definition of the spatial reference system.
• gdalinfo on the file still reports "Coordinate System is `'"
• opening the file with the python bindings and calling dst.GetProjection() returns a blank string (i.e. no projection information)

The full listings of my spatial reference definition, the output of ncdump and gdalinfo and the output from the python session are listed below. Can anyone tell me whether I'm doing something wrong, or should I not count on georeferencing being supported by the NetCDF driver?

An example NetCDF file with the global attribute "spatial_ref" is attached.

~James

---

james@basalt ~ $ ncatted -a spatial_ref,global,a,c,'PROJCS["unnamed",GEOGCS["Normal Sphere
(r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]'
narccap_ecpc-20c3m-tasmin-ncep.nc narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc

---

james@basalt ~ $ ncdump -h narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
netcdf narccap_ecpc-20c3m-tasmin-ncep-with-srs {
dimensions:
	xc = 123 ;
	yc = 104 ;
	time = UNLIMITED ; // (9496 currently)
	bnds = 2 ;
variables:
	char polar_stereographic ;
		polar_stereographic:grid_mapping_name =
"polar_stereographic" ;
		polar_stereographic:longitude_of_central_meridian = 263. ;
		polar_stereographic:straight_vertical_longitude_from_pole =
263. ;
		polar_stereographic:standard_parallel = 60. ;
		polar_stereographic:TrueGridLength_Latitude = 60. ;
		polar_stereographic:XGridLength_M = 50000. ;
		polar_stereographic:YGridLength_M = 50000. ;
		polar_stereographic:PoleOnPlane = "north" ;
		polar_stereographic:false_easting = 3925000. ;
		polar_stereographic:false_northing = 8325000. ;
		polar_stereographic:latitude_of_projection_origin = 90. ;
	double yc(yc) ;
		yc:long_name = "y-coordinate of Cartesian system" ;
		yc:standard_name = "projection_y_coordinate" ;
		yc:axis = "Y" ;
		yc:units = "m" ;
	double xc(xc) ;
		xc:long_name = "x-coordinate of Cartesian system" ;
		xc:standard_name = "projection_x_coordinate" ;
		xc:axis = "X" ;
		xc:units = "m" ;
	double lon(yc, xc) ;
		lon:units = "degrees_east" ;
		lon:long_name = "longitude" ;
		lon:standard_name = "longitude" ;
		lon:axis = "X" ;
		lon:NX_NumPntsOnParallel = 123 ;
		lon:NY_NumPntsOnMeridian = 104 ;
		lon:actual_range = 211.5412f, 316.2759f ;
		lon:XGridLength_M = 50000. ;
		lon:YGridLength_M = 50000. ;
		lon:PoleOnPlane = "north" ;
	double lat(yc, xc) ;
		lat:units = "degrees_north" ;
		lat:long_name = "latitude" ;
		lat:standard_name = "latitude" ;
		lat:axis = "Y" ;
		lat:NX_NumPntsOnParallel = 123 ;
		lat:NY_NumPntsOnMeridian = 104 ;
		lat:actual_range = 21.23807f, 67.63758f ;
		lat:XGridLength_M = 50000. ;
		lat:YGridLength_M = 50000. ;
		lat:PoleOnPlane = "north" ;
	double level ;
		level:long_name = "height" ;
		level:units = "m" ;
		level:positive = "up" ;
	double time(time) ;
		time:long_name = "time" ;
		time:standard_name = "time" ;
		time:axis = "T" ;
		time:units = "days since 1979-01-01 00:00:0.0" ;
		time:calendar = "gregorian" ;
		time:bounds = "time_bnds" ;
	double time_bnds(time, bnds) ;
	float tasmin(time, yc, xc) ;
		tasmin:units = "K" ;
		tasmin:standard_name = "air_temperature" ;
		tasmin:long_name = "Minimum Daily Surface Air Temperature" ;
		tasmin:cell_methods = "time: minimum (interval: 1 day)" ;
		tasmin:missing_value = 1.e+20f ;
		tasmin:coordinates = "lon lat level" ;
		tasmin:grid_mapping = "polar_stereographic" ;
		tasmin:grid_mapping_name = "polar_stereographic" ;
		tasmin:type = "Instantaneous" ;
		tasmin:actual_min = 206.81f ;
		tasmin:actual_max = 307.07f ;
		tasmin:NumMissing = 0 ;
		tasmin:_FillValue = 1.e+20f ;
		tasmin:NX_NumPntsOnParallel = 123 ;
		tasmin:NY_NumPntsOnMeridian = 104 ;
		tasmin:LowerLeftLatitude = 21.238 ;
		tasmin:LowerLeftLongitude = 211.541 ;
		tasmin:XGridLength_M = 50000. ;
		tasmin:YGridLength_M = 50000. ;
		tasmin:PoleOnPlane = "north" ;



// global attributes:
		:Conventions = "CF-1.0" ;
		:institution = "Experimental Climate Prediction Center
(ECPC), Scripps Institution of Oceanography, University of California, San
Diego, La Jolla, California, USA" ;
		:source = "ECPC Regional Spectral Model
(2005):atmosphere:ECPC-RSM hydrostatic; 28-sigma vertical layers; land:Noah"
;
		:project_id = "NARCCAP" ;
		:table_id = "Table 1" ;
		:experiment_id = "present-day climate using NCEP/DOE
		Reanalysis" ;
		:title = "ECPC-RSM output prepared for NARCCAP present-day
		climate using NCEP/DOE Reanalyis" ;
		:realization = "1" ;
		:contact1 = "A. Nunes" ;
		:contact2 = "anunes@ucsd.edu" ;
		:history1 = "Original output in GRIB 1" ;
		:history2 = "Contains modifications made to original data" ;
		:history = "Mon Aug  9 13:44:50 2010: ncatted -a
		spatial_ref,global,a,c,PROJCS[\"unnamed\",GEOGCS[\"Normal
		Sphere
		(r=6370997)\",DATUM[\"unknown\",SPHEROID[\"sphere\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",263],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",3925000],PARAMETER[\"false_northing\",8325000]]
		narccap_ecpc-20c3m-tasmin-ncep.nc
		narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc\n",
			"Fri May  8 09:34:14 2009: ncrcat
tasmin_ECPC_1979010106.nc tasmin_ECPC_1981010106.nc
tasmin_ECPC_1986010106.nc tasmin_ECPC_1991010106.nc
tasmin_ECPC_1996010106.nc tasmin_ECPC_2001010106.nc
narccap_ecpc-20c3m-tasmin-ncep.nc" ;
		:nco_openmp_thread_number = 1 ;
		:spatial_ref = "PROJCS[\"unnamed\",GEOGCS[\"Normal Sphere
(r=6370997)\",DATUM[\"unknown\",SPHEROID[\"sphere\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",263],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",3925000],PARAMETER[\"false_northing\",8325000]]"
;
}

---

james@basalt ~ $ gdalinfo narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
Driver: netCDF/Network Common Data Format
Files: narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
       narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc.aux.xml
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#institution=Experimental Climate Prediction Center (ECPC),
Scripps Institution of Oceanography, University of California, San Diego, La
Jolla, California, USA
  NC_GLOBAL#source=ECPC Regional Spectral Model (2005):atmosphere:ECPC-RSM
hydrostatic; 28-sigma vertical layers; land:Noah
  NC_GLOBAL#project_id=NARCCAP
  NC_GLOBAL#table_id=Table 1
  NC_GLOBAL#experiment_id=present-day climate using NCEP/DOE Reanalysis
  NC_GLOBAL#title=ECPC-RSM output prepared for NARCCAP present-day climate
  using NCEP/DOE Reanalyis
  NC_GLOBAL#realization=1
  NC_GLOBAL#contact1=A. Nunes
  NC_GLOBAL#contact2=anunes@ucsd.edu
  NC_GLOBAL#history1=Original output in GRIB 1
  NC_GLOBAL#history2=Contains modifications made to original data
  NC_GLOBAL#history=Mon Aug  9 13:44:50 2010: ncatted -a
  spatial_ref,global,a,c,PROJCS["unnamed",GEOGCS["Normal Sphere
  (r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]
  narccap_ecpc-20c3m-tasmin-ncep.nc
  narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
Fri May  8 09:34:14 2009: ncrcat tasmin_ECPC_1979010106.nc
tasmin_ECPC_1981010106.nc tasmin_ECPC_1986010106.nc
tasmin_ECPC_1991010106.nc tasmin_ECPC_1996010106.nc
tasmin_ECPC_2001010106.nc narccap_ecpc-20c3m-tasmin-ncep.nc
  NC_GLOBAL#nco_openmp_thread_number=1
  NC_GLOBAL#spatial_ref=PROJCS["unnamed",GEOGCS["Normal Sphere
(r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc":lon
  SUBDATASET_1_DESC=[104x123] longitude (64-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc":lat
  SUBDATASET_2_DESC=[104x123] latitude (64-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc":time_bnds
  SUBDATASET_3_DESC=[9496x2] time_bnds (64-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc":tasmin
  SUBDATASET_4_DESC=[9496x104x123] air_temperature (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

---

>>> from osgeo import gdal
>>> dst = gdal.Open("narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc")
>>> dst.GetProjection()
''

[hiebert]

Example netcdf which has a global attribute "spatial_ref" containing WKT that describes the spatial reference system

[Kyle Shannon]

I added support for polar stereographic in r 18931, as well as 3 other projections. These were not back ported to 1.7. I will test this dataset and post results tomorrow.

[Kyle Shannon]

Frank, I will take ownership of this, as I think I already added support. We may want to backport a few of my commits to add the support for this ticket and #3520. I have to look it over and make sure I commited correctly.

[ilucena]

Hi guys,

It looks like it just a matter of informing the subdataset properly. See:

gdalinfo NETCDF:narccap_ecpc-20c3m-tasmin-ncep-with-srs-one-t.nc:tasmin
Driver: netCDF/Network Common Data Format
Files: narccap_ecpc-20c3m-tasmin-ncep-with-srs-one-t.nc
Size is 123, 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["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",3925000],
    PARAMETER["false_northing",8325000]]
Origin = (950000.000000000000000,6000000.000000000000000)
Pixel Size = (50000.000000000000000,-50000.000000000000000)
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#institution=Experimental Climate Prediction Center (ECPC), Scripps Institution of Oceanography, University of California, San Diego, La Jolla, California, USA
  NC_GLOBAL#source=ECPC Regional Spectral Model (2005):atmosphere:ECPC-RSM hydrostatic; 28-sigma vertical layers; land:Noah
  NC_GLOBAL#project_id=NARCCAP
  NC_GLOBAL#table_id=Table 1
  NC_GLOBAL#experiment_id=present-day climate using NCEP/DOE Reanalysis
  NC_GLOBAL#title=ECPC-RSM output prepared for NARCCAP present-day climate using NCEP/DOE Reanalyis
  NC_GLOBAL#realization=1
  NC_GLOBAL#contact1=A. Nunes
  NC_GLOBAL#contact2=anunes@ucsd.edu
  NC_GLOBAL#history1=Original output in GRIB 1
  NC_GLOBAL#history2=Contains modifications made to original data
  NC_GLOBAL#history=Mon Aug  9 16:38:25 2010: ncks -d time,1,1 narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc narccap_ecpc-20c3m-tasmin-ncep-with-srs-one-t.nc
Mon Aug  9 13:44:50 2010: ncatted -a spatial_ref,global,a,c,PROJCS["unnamed",GEOGCS["Normal Sphere (r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]] narccap_ecpc-20c3m-tasmin-ncep.nc narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
Fri May  8 09:34:14 2009: ncrcat tasmin_ECPC_1979010106.nc tasmin_ECPC_1981010106.nc tasmin_ECPC_1986010106.nc tasmin_ECPC_1991010106.nc tasmin_ECPC_1996010106.nc tasmin_ECPC_2001010106.nc narccap_ecpc-20c3m-tasmin-ncep.nc
  NC_GLOBAL#nco_openmp_thread_number=1
  NC_GLOBAL#spatial_ref=PROJCS["unnamed",GEOGCS["Normal Sphere (r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]
  NC_GLOBAL#NCO=3.9.9
  tasmin#units=K
  tasmin#standard_name=air_temperature
  tasmin#long_name=Minimum Daily Surface Air Temperature
  tasmin#cell_methods=time: minimum (interval: 1 day)
  tasmin#missing_value=1.000000e+20
  tasmin#coordinates=lon lat level
  tasmin#grid_mapping=polar_stereographic
  tasmin#grid_mapping_name=polar_stereographic
  tasmin#type=Instantaneous
  tasmin#actual_min=2.068100e+02
  tasmin#actual_max=3.070700e+02
  tasmin#NumMissing=0
  tasmin#_FillValue=1.000000e+20
  tasmin#NX_NumPntsOnParallel=123
  tasmin#NY_NumPntsOnMeridian=104
  tasmin#LowerLeftLatitude=21.238
  tasmin#LowerLeftLongitude=211.541
  tasmin#XGridLength_M=50000
  tasmin#YGridLength_M=50000
  tasmin#PoleOnPlane=north
  yc#long_name=y-coordinate of Cartesian system
  yc#standard_name=projection_y_coordinate
  yc#axis=Y
  yc#units=m
  xc#long_name=x-coordinate of Cartesian system
  xc#standard_name=projection_x_coordinate
  xc#axis=X
  xc#units=m
  time#long_name=time
  time#standard_name=time
  time#axis=T
  time#units=days since 1979-01-01 00:00:0.0
  time#calendar=gregorian
  time#bounds=time_bnds
  polar_stereographic#grid_mapping_name=polar_stereographic
  polar_stereographic#longitude_of_central_meridian=263
  polar_stereographic#straight_vertical_longitude_from_pole=263
  polar_stereographic#standard_parallel=60
  polar_stereographic#TrueGridLength_Latitude=60
  polar_stereographic#XGridLength_M=50000
  polar_stereographic#YGridLength_M=50000
  polar_stereographic#PoleOnPlane=north
  polar_stereographic#false_easting=3.925e+06
  polar_stereographic#false_northing=8.325e+06
  polar_stereographic#latitude_of_projection_origin=90
Corner Coordinates:
Upper Left  (  950000.000, 6000000.000) ( 51d59'30.93"W, 57d 5'39.85"N)
Lower Left  (  950000.000,  800000.000) ( 21d34'16.71"W, 25d11'50.63"N)
Upper Right ( 7100000.000, 6000000.000) ( 53d47'7.08"E, 55d46'52.23"N)
Lower Right ( 7100000.000,  800000.000) ( 22d52'34.28"E, 24d42'40.77"N)
Center      ( 4025000.000, 3400000.000) (  1d 9'47.54"E, 47d49'57.67"N)
Band 1 Block=123x1 Type=Float32, ColorInterp=Undefined
  NoData Value=1.00000002004087734e+20
  Metadata:
    NETCDF_VARNAME=tasmin
    NETCDF_DIMENSION_time=1.25
    NETCDF_time_units=days since 1979-01-01 00:00:0.0

[Kyle Shannon]

That makes sense, although I would like to take a look into the polar stereographic portion of this. I am pretty sure there was no explicit setting of polar stereographic before r18931. Thanks for the bug report.

[Kyle Shannon]

Are you sure of your version? I just built 1.7.2 and received a for the SRS. Can you post what gdalinfo --version reports? It worked fine on my 1.8dev build.

[Kyle Shannon]

After looking at the file a bit, I realized that the spatial_ref tag was in the test file. gdal was reading this in as wkt and setting the projection. 1.8dev reads the cf-1.x grid_mapping tags and sets the projection. If you would like to read the projection information from the file without adding the gdal tags, you will have to update to the trunk. I am closing this ticket, re-open it if you see fit.

[hiebert]

So just to make sure that I understand correctly:

The reason that I don't get an SRS in 1.7.2 is that even though GDAL is reading the spatial_ref tag, the polar stereographic project wasn't supported? But I should be able to expect all of the other projections listed in static const oNetcdfSRS poNetcdfSRS in netcdfdataset.h to work?

And if I use 1.8dev, then I shouldn't even need to add the spatial_ref tag, because GDAL will properly read the NetCDF metadata from the cf grid_mapping tags?

Many thank for taking a look at this. It was very helpful.

[Kyle Shannon]

Sorry it took me so long to get back. So,

""The reason that I don't get an SRS in 1.7.2 is that even though GDAL is reading the spatial_ref tag, the polar stereographic project wasn't supported? ""

No, that isn't correct. If you take a peek at ilucena's comment above about the subdatasets, In order to get the georeferencing, you have to call one of the subdatasets.

""But I should be able to expect all of the other projections listed in static const oNetcdfSRS poNetcdfSRS in netcdfdataset.h to work?""

In reference to the spatial_ref tag, I do believe they work. but I am not sure, that was before my time, I will have to look into it.

""And if I use 1.8dev, then I shouldn't even need to add the spatial_ref tag, because GDAL will properly read the NetCDF metadata from the cf grid_mapping tags?""

Correct. GDAL should now support grid_mapping tags in cf 1.4. I haven't checked to see if they added any new grid mappings in 1.5. I will do that today.

Thanks. Go ahead and reopen this if it is unclear, or post a question again.

kss

[hiebert]

Replying to kyle:

Sorry it took me so long to get back. So,

No problem.

No, that isn't correct. If you take a peek at ilucena's comment above about the subdatasets, In order to get the georeferencing, you have to call one of the subdatasets.

Hmmm, I have tried it with what I'm pretty sure is the correct subdataset in 1.7.2 and I don't get a coordinate system.

$ gdalinfo --version
GDAL 1.7.2, released 2010/04/23
$ gdalinfo NETCDF:narccap_ecpc-20c3m-tasmin-ncep-with-srs-one-t.nc:tasmin | grep Coord
Coordinate System is `'

That said, don't burn too much (or any) time on this. We have a reasonable work-around for our application at present, and as long as we'll have direct support of the grid_mapping tags in 1.8.0, I'm happy. I just tried it with rev21059 from 1.8dev and I get a CRS, so it should be all good.

Correct. GDAL should now support grid_mapping tags in cf 1.4. I haven't checked to see if they added any new grid mappings in 1.5. I will do that today.

Fantastic. I appreciate the help!