Opened 14 years ago

Closed 10 years ago

Last modified 10 years ago

#1637 closed enhancement (fixed)

patch for 1.4.x : frmts/netcdf (NC_GLOBAL#{spatial_ref/GeoTransform})

Reported by: cermak@… Owned by: dnadeau
Priority: normal Milestone: 1.8.1
Component: GDAL_Raster Version: 1.4.1
Severity: normal Keywords: netcdf
Cc: cermak@…, warmerdam, Mateusz Łoskot, Markus Neteler, Kyle Shannon

Description (last modified by Mateusz Łoskot)

Recent events have us looking the the direct access Netcdf driver (frmts/netcdf). We discovered that the spatial_ref and GeoTransform? global tags are ignored for non-complient CF-1.0 files.

Patch: http://ak.aoos.org/data/patch/frmts_netcdf_global_georef.patch

WRF data file, if needed: http://137.229.130.18/data/forecasts/AEFF/WRF/CIPWS12km/200705/14/00/wrfout_d01_2007-05-16_00.nc

(sorry about the size! 14 MB)

1.4.1:

gdalinfo NETCDF:"wrfout_d01_2007-05-16_00.nc":T2
  NC_GLOBAL#spatial_ref=+proj=stere +lon_0=-147.300003052 +lat_0=90
+lat_ts=60.2000007629 +a=6378273 +rf=298.279405043 +x_0=0 +y_0=0
  NC_GLOBAL#GeoTransform=-541931.561016 12045.5692702 0 -3610575.06283 0
12025.6894815
  T2#FieldType=104
  T2#MemoryOrder=XY
  T2#description=TEMP at 2 M
  T2#units=K
  T2#stagger=
  T2#coordinates=XLONG XLAT
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,   72.0)
Upper Right (   90.0,    0.0)
Lower Right (   90.0,   72.0)
Center      (   45.0,   36.0)
Band 1 Block=90x1 Type=Float32, ColorInterp=Undefined
  NoData Value=9.96920996838687e+36
  Metadata:
    NETCDF_VARNAME=T2
    NETCDF_DIMENSION_Time=1
GDAL: GDALClose(NETCDF:wrfout_d01_2007-05-16_00.nc:T2)
GDAL: GDALDeregister_GTiff() called.

With the above patch:

  NC_GLOBAL#spatial_ref=+proj=stere +lon_0=-147.300003052 +lat_0=90
+lat_ts=60.2000007629 +a=6378273 +rf=298.279405043 +x_0=0 +y_0=0
  NC_GLOBAL#GeoTransform=-541931.561016 12045.5692702 0 -3610575.06283 0
12025.6894815
  T2#FieldType=104
  T2#MemoryOrder=XY
  T2#description=TEMP at 2 M
  T2#units=K
  T2#stagger=
  T2#coordinates=XLONG XLAT
Corner Coordinates:
Upper Left  ( -541931.561,-3610575.063)
Lower Left  ( -541931.561,-2744725.420)
Upper Right (  542169.673,-3610575.063)
Lower Right (  542169.673,-2744725.420)
Center      (     119.056,-3177650.241)
Band 1 Block=90x1 Type=Float32, ColorInterp=Undefined
  NoData Value=9.96920996838687e+36
  Metadata:
    NETCDF_VARNAME=T2
    NETCDF_DIMENSION_Time=1
GDAL: GDALClose(NETCDF:wrfout_d01_2007-05-16_00.nc:T2)
GDAL: GDALDeregister_GTiff() called.

Change History (12)

comment:1 Changed 14 years ago by guest

Component: defaultGDAL_Raster
Owner: changed from warmerdam to cermak@…

Unfortunately, this only seems 1/2 the battle. I was hoping the dataset would then be accessible via Mapserver as we are able to do via OPeNDAP

I will reassign to myself until I can figure out the 2nd half.

I copied a bit of code to convert the spatial_ref if it was a PROJ.4 string to WKT. gdalwarp didn't like the PROJ.4 string we use.

Version 2 vs. 1.4.1: http://ak.aoos.org/data/patch/frmts_netcdf_global_georef.patch

To crack the 2nd half of the nut, I need to get a working example from the OPeNDAP side. Right now, I get a blank image.

gdalwarp NETCDF:/space/data/forecasts/AEFF/WRF/CIPWS12km/200705/14/00/wrfout_d01_2007-05-16_00.nc:T2 x.tiff GDAL: GDALOpen(x.tiff) succeeds as GTiff. GDAL_netCDF: dim_count = 9

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute GDAL_netCDF: var_count = 130

GDAL_netCDF: Found NC_GLOBAL#spatial_ref

GDAL_netCDF: Fetching NC_GLOBAL#spatial_ref

GDAL_netCDF: WKT = PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6378273,298.279405043]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTIONPolar_Stereographic?,PARAMETER["latitude_of_origin",60.2000007629],PARAMETER["central_meridian",-147.300003052],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0]]

GDAL_netCDF: Fetching NC_GLOBAL#GeoTransform?

GDAL: GDALOpen(NETCDF:/space/data/forecasts/AEFF/WRF/CIPWS12km/200705/14/00/wrfout_d01_2007-05-16_00.nc:T2) succeeds as netCDF. Processing input file NETCDF:/space/data/forecasts/AEFF/WRF/CIPWS12km/200705/14/00/wrfout_d01_2007-05-16_00.nc:T2. GDAL: GDALWarpKernel()::GWKNearestNoMasksFloat() Src=0,0,90x72 Dst=0,0,90x72 :0...11...20...30...40...50...61...70...80...90...100 - done. GDAL: GDALClose(NETCDF:/space/data/forecasts/AEFF/WRF/CIPWS12km/200705/14/00/wrfout_d01_2007-05-16_00.nc:T2) GDAL: GDALClose(x.tiff) GDAL: GDALDeregister_GTiff() called.

comment:2 Changed 14 years ago by guest

Owner: changed from cermak@… to warmerdam

Ok. As it turns out, the patch works fine.

I happen to have a THREDDS OPeDNAP server to help me reproduce the graphics for both drivers.

The one major difference is that the Netcdf driver does not allow selection of bands like the OPeNDAP driver. Typically, we have all the time dimensions stacked together so we can do something like: T2?[3][y][x] where 3 points to the timeslice we want out of the dataset. If there are multiple times, the NETCDF driver likes to create multiple bands and then the behavior gets a little wierd... the WKT gets stomped somehow.

Creating output file that is 319P x 224L. GDAL: GDALDriver::Create(GTiff,test.tif,319,224,17,Int16,(nil)) ERROR 1: Only OGC WKT Projections supported for writing to GeoTIFF. (X17�unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6378273,298.279405043]],P RIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTIONPolar_Stereographic?,PARAMETER ["latitude_of_origin",67.5],PARAMETER["central_meridian",-150],PARAMETER["scale_factor",1],PARAMETE R["false_easting",0],PARAMETER["false_northing",0]] not supported.

Overall, this patch will allow georeferencing to work if the global attributes are defined. I think someone else noticed that this stopped working as a couple versions ago.

comment:3 Changed 14 years ago by warmerdam

Cc: warmerdam added
Milestone: 1.5.0
Owner: changed from warmerdam to dnadeau
Reporter: changed from guest to cermak@…

Denis,

Can you review this and consider applying the patch? I'm a bit lost.

comment:4 in reply to:  3 Changed 14 years ago by guest

Replying to warmerdam:

Denis,

Can you review this and consider applying the patch? I'm a bit lost.

The patch allows spatial_ref and GeoTransform? work via NC_GLOBAL metadata tags. As of 1.4.1, even with the metadata tags, it would not georeference properly.

comment:5 Changed 14 years ago by Mateusz Łoskot

Cc: Mateusz Łoskot added

comment:6 Changed 13 years ago by Mateusz Łoskot

Description: modified (diff)
Keywords: netcdf added

comment:7 Changed 13 years ago by Markus Neteler

Cc: Markus Neteler added

I have applied the patch (having the same problem of no spatial_ref with netCDF files from NSIDC (eg, ftp://sidads.colorado.edu/pub/DATASETS/snow/nsidc0321v01/north/NL.2001.nsidc0321v01.nc). The files are coming in EASE-Grid but it is not detected by GDAL, even with this patch:

gdalinfo /geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: /geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#creation_date=Wed Feb 20 13:14:48 MST 2008
  NC_GLOBAL#grid_id=Nl
  NC_GLOBAL#references=Data set documentation: http://nsidc.org/data/nsidc-0321.html
  NC_GLOBAL#comment=Not set at this time
  NC_GLOBAL#source=Satellite observations from DMSP SSM/I and Terra MODIS
  NC_GLOBAL#institution=National Snow & Ice Data Center, Boulder, CO
  NC_GLOBAL#title=Global EASE-Grid 8-day Blended SSM/I and MODIS Snow Cover
  NC_GLOBAL#history=Wed Feb 20 13:15:33 2008: ncrcat -O NL.2001001.nsidc0321v01.nc NL.2001009.nsidc0321v01.nc NL.2001017.nsidc0321v01.nc NL.2001025.nsidc0321v01.nc NL.2001033.nsidc0321v01.nc NL.2001041.nsidc0321v01.nc NL.2001049.nsidc0321v01.nc NL.2001057.nsidc0321v01.nc NL.2001065.nsidc0321v01.nc NL.2001073.nsidc0321v01.nc NL.2001081.nsidc0321v01.nc NL.2001089.nsidc0321v01.nc NL.2001097.nsidc0321v01.nc NL.2001105.nsidc0321v01.nc NL.2001113.nsidc0321v01.nc NL.2001121.nsidc0321v01.nc NL.2001129.nsidc0321v01.nc NL.2001137.nsidc0321v01.nc NL.2001145.nsidc0321v01.nc NL.2001153.nsidc0321v01.nc NL.2001161.nsidc0321v01.nc NL.2001185.nsidc0321v01.nc NL.2001193.nsidc0321v01.nc NL.2001201.nsidc0321v01.nc NL.2001217.nsidc0321v01.nc NL.2001225.nsidc0321v01.nc NL.2001233.nsidc0321v01.nc NL.2001241.nsidc0321v01.nc NL.2001249.nsidc0321v01.nc NL.2001257.nsidc0321v01.nc NL.2001265.nsidc0321v01.nc NL.2001273.nsidc0321v01.nc NL.2001281.nsidc0321v01.nc NL.2001289.nsidc0321v01.nc NL.2001297.nsidc0321v01.nc NL.2001305.nsidc0321v01.nc NL.2001313.nsidc0321v01.nc NL.2001321.nsidc0321v01.nc NL.2001329.nsidc0321v01.nc NL.2001337.nsidc0321v01.nc NL.2001345.nsidc0321v01.nc NL.2001353.nsidc0321v01.nc NL.2001361.nsidc0321v01.nc NL.2001.nsidc0321v01.nc
  NC_GLOBAL#nco_openmp_thread_number=1
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"/geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc":bpInfo
  SUBDATASET_1_DESC=[43x81] bpInfo ()
  SUBDATASET_2_NAME=NETCDF:"/geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc":scaInfo
  SUBDATASET_2_DESC=[43x81] scaInfo ()
  SUBDATASET_3_NAME=NETCDF:"/geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc":latitude
  SUBDATASET_3_DESC=[721x721] latitude (64-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"/geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc":longitude
  SUBDATASET_4_DESC=[721x721] longitude (64-bit floating-point)
  SUBDATASET_5_NAME=NETCDF:"/geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc":SWE
  SUBDATASET_5_DESC=[43x721x721] lwe_thickness_of_surface_snow_amount (8-bit integer)
  SUBDATASET_6_NAME=NETCDF:"/geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc":SCA
  SUBDATASET_6_DESC=[43x721x721] surface_snow_area_fraction (8-bit integer)
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)
Warning 1: Unable to save auxilary information in /geodata2_originals_raid5/world_snowcover_2000_today_weekly_0_28deg/north/NL.2001.nsidc0321v01.nc.aux.xml.

Note - this may need a fix:

netcdfdataset.cpp:1368: warning: 'nVarID' may be used uninitialized in this function

Best, Markus

comment:8 Changed 12 years ago by Kyle Shannon

Cc: Kyle Shannon added

comment:9 Changed 12 years ago by Kyle Shannon

Are you writing the tags during post-processing of wrf?

comment:10 Changed 12 years ago by guest

No. At the time of this ticket, we ran the WRF output through another post-processing script (python) to add the necessary attribute tags to enable gdal to produce the proper georeferencing information.

comment:11 Changed 10 years ago by etourigny

Resolution: fixed
Status: newclosed

Example file (WRF) and patch are no longer accessible.

The recent driver (in trunk) reads the EASE-Grid files from NSIDC (in LAEA projection), using the CF atributes:

 $ gdalinfo NETCDF:NL.2001.nsidc0321v01.nc:SWE
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: none associated
Size is 721, 721
Coordinate System is:
PROJCS["LAEA (WGS84) ",
    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["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",90],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-0.500000000000000,720.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
....

ncdump output (with CF projection variable)

$ ncdump -h  NL.2001.nsidc0321v01.nc 
netcdf NL.2001.nsidc0321v01 {
dimensions:
	time = UNLIMITED ; // (43 currently)
	rows = 721 ;
	cols = 721 ;
	maxFilenameLength = 81 ;
	ncl4 = 10 ;
variables:
	double time(time) ;
		time:units = "days since -4713-01-01 12:00:00" ;
		time:long_name = "julian day" ;
...
	short SWE(time, rows, cols) ;
		SWE:standard_name = "lwe_thickness_of_surface_snow_amount" ;
		SWE:grid_mapping = "projection" ;
		SWE:ancillary_variables = "SCA bpInfo scaInfo" ;
		SWE:flag_meanings = "\n",
			"no_microwave_data \n",
			"grid_corners \n",
			"ocean \n",
			"permanent_ice_sheets_and_large_glaciers\n",
			"no_microwave_SWE_but_visible_SCA_gt_25_pcent" ;
		SWE:flag_values = -150s, -200s, -250s, -300s, -350s ;
		SWE:_FillValue = -200s ;
		SWE:information = "\n",
			" > 0        : SWE from deep microwave algorithm for 8-day period\n",
			" 0          : no snow\n",
			" -100 to -1 : -1 * SWE from shallow microwave algorithm\n",
			" -150       : missing microwave brightness temperatures\n",
			" -200       : fixed value for corners\n",
			" -250       : ocean\n",
			" -300       : permanent ice sheets and large glaciers\n",
			" -350       : no microwave SWE, but visible SCA > 25%" ;
		SWE:coordinates = "latitude longitude" ;
		SWE:units = "mm" ;
		SWE:long_name = "Snow Water Equivalent" ;
...
	char projection(ncl4) ;
		projection:grid_id = "NL" ;
		projection:reference = "EASE-Grid projection definition documention: http://nsidc.org/data/ease/ease_grid.html" ;
		projection:NSIDC_mapx_grid_parameter_definition_file = "Nl.gpd" ;
		projection:grid_mapping_name = "lambert_azimuthal_equal_area" ;
		projection:longitude_of_projection_origin = 0.f ;
		projection:latitude_of_projection_origin = 90.f ;
		projection:false_easting = 0.f ;
		projection:false_northing = 0.f ;
		projection:projection_x_coordinate_origin = 360.f ;
		projection:projection_y_coordinate_origin = 360.f ;
		projection:scale_factor_at_projection_origin = 25.06753f ;

// global attributes:
		:creation_date = "Wed Feb 20 13:14:48 MST 2008" ;
...
}

Closing this bug as it is outdated and trunk can read the NSIDC files.

comment:12 Changed 10 years ago by etourigny

as reported in bug #2584 the netcdf files from NSIDC containt row/column numbers instead of LAEA coordinates. The current netcdf driver is not able to deal with these values.

Note: See TracTickets for help on using tickets.