James Hiebert's Blog

GDAL hates netCDF

Unfortunately I can’t seem to figure out how to get GDAL to read in a netcdf file, and detect the spatial reference system.

This page (in the section titled “Georeference”) claims that if the file has an attribute named “spatial_ref”, it can read it and set up the georeferencing. OK, let’s try it.

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]]"
;
}

OK, it has the attribute. Let’s take a look at what gdal says about it:

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)

Note that the Coordinate System is still blank. Damn, what’s going on?

Update (2010-08-20):

Turns out that GDAL doesn’t hate netcdf. After working with the gdal-dev mailing list and submitting a ticket to OSgeo, we figured out that GDAL was reading the metadata correctly, but polar stereographic projection just hadn’t yet been supported in the NetCDF driver. And even better news was that in future versions of GDAL, CF-metadata tags (i.e. the tags on all of our model data) will be supported directly. Great news!

blog comments powered by Disqus

Published

09 August 2010

Category

work

Tags