[gdal-dev] How does GDAL's NetCDF driver determine the spatial
reference system
James Hiebert
hiebert at uvic.ca
Mon Aug 9 18:21:51 EDT 2010
Hi All,
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?
~James
james at 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 at 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" ;
...
// global attributes:
...
: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 at 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#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)
james at basalt ~ $ python
>>> from osgeo import gdal
>>> dst = gdal.Open("narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc")
>>> dst.GetProjection()
More information about the gdal-dev
mailing list