Opened 11 years ago

Closed 10 years ago

#3863 closed defect (invalid)

Error while reprojecting with gdalwarp an ECMWF netCDF

Reported by: toyzerocha Owned by: Kyle Shannon
Priority: normal Milestone: 1.8.1
Component: GDAL_Raster Version: 1.7.1
Severity: normal Keywords: netcdf, projection
Cc: Kyle Shannon

Description (last modified by Kyle Shannon)

I'm trying to reproject an ECMWF netCDF file to a specific coordinate system EPSG:3726 with the following command:

gdalwarp -t_srs EPSG:3726 2009.nc 2009.tif

and I projects to this coordinate system but the coordinates are not changed. This means that, the coordinates are copied from the original system to this one without any reprojection.

This happened in GDAL Version 1.7.1 in Ubuntu and Windows

This is related with the following messages: http://lists.osgeo.org/pipermail/gdal-dev/2010-December/026956.html http://lists.osgeo.org/pipermail/gdal-dev/2010-December/026954.html

Attachments (2)

2009.nc (20.7 KB) - added by toyzerocha 11 years ago.
NetCDF original file
2009NCwindows.tif (21.8 KB) - added by toyzerocha 11 years ago.
Tif file after gdalwarp

Download all attachments as: .zip

Change History (7)

Changed 11 years ago by toyzerocha

Attachment: 2009.nc added

NetCDF original file

Changed 11 years ago by toyzerocha

Attachment: 2009NCwindows.tif added

Tif file after gdalwarp

comment:1 Changed 11 years ago by Even Rouault

Cc: Kyle Shannon added

comment:2 Changed 11 years ago by Kyle Shannon

Description: modified (diff)

comment:3 Changed 11 years ago by Kyle Shannon

Description: modified (diff)
Milestone: 1.8.0
Owner: changed from warmerdam to Kyle Shannon
Severity: majornormal
Status: newassigned

I apologize ahead of time for the long message. The short answer is:

  • The netcdf driver does not have enough information to specify a spatial reference
  • The data is obviously geographic, but on what sphere with what prime meridian?
  • Even if WGS84 is assumed, the warp fails due to a failure in proj4/gdal to warp EPSG:3726 to EPSG:4326
  • Partial support for a partial solution is in 1.8dev, but definitely not all.
  • See below for a little more info

I have been looking at bits of this issue. There are a few things that jump out at me. I notice that the netcdf file contains latitude and longitude as variables. I do believe the driver attempts to geo-reference this, but there are some issues. The driver picks up the geotransform values, meaning origin and pixel size:

gdalinfo 2009.nc -nomd 
Driver: netCDF/Network Common Data Format
Files: 2009.nc
Size is 3, 4
Coordinate System is `'
Origin = (350.250000000000000,42.750000000000000)
Pixel Size = (1.500000000000000,-1.500000000000000)
Corner Coordinates:
Upper Left  (     350.250,      42.750) 
Lower Left  (     350.250,      36.750) 
Upper Right (     354.750,      42.750) 
Lower Right (     354.750,      36.750) 
Center      (     352.500,      39.750) 

but obviously no spatial reference. This is because the driver cannot specify the parameters of the sphere or datum. It knows which latitude and longitude to use, but it doesn't know on what sphere or what meridian to reference.

If I add the following lines to the netcdf file using ncdump/emacs/ncgen:

tp:long_name = "Total precipitation" ;
tp:grid_mapping = "lat_lon_map" ;
char lat_lon_map ;
	     	lat_lon_map:grid_mapping_name = "latitude_longitude" ;

The driver reports(1.8dev/trunk):

gdalinfo new.nc -nomd
Coordinate System is:
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"]]
Origin = (350.250000000000000,42.750000000000000)
Pixel Size = (1.500000000000000,-1.500000000000000)
Corner Coordinates:
Upper Left  (     350.250,      42.750) (350d15' 0.00"E, 42d45' 0.00"N)
Lower Left  (     350.250,      36.750) (350d15' 0.00"E, 36d45' 0.00"N)
Upper Right (     354.750,      42.750) (354d45' 0.00"E, 42d45' 0.00"N)
Lower Right (     354.750,      36.750) (354d45' 0.00"E, 36d45' 0.00"N)
Center      (     352.500,      39.750) (352d30' 0.00"E, 39d45' 0.00"N)

Even this is suspect, as it defaults to well known geographic system, WGS84. It probably shouldn't do that. The support for the CF-1.x tags used above is available in 1.8dev, not 1.7.x. According to CF-1.x, if the latitude_longitude grid_mapping is defined, the sphere should be defined. It isn't so the driver sets it to WGS-84.

Even after this assigning a generic WGS84 reference, gdalwarp fails. It appears gdal/proj4 has trouble converting from ESPG:3726 to WGS84. I am not exactly sure where it needs to do this, but I imagine it is to check bounds of the output image. It also looks like UTM 19 is quite a ways away from ~350 deg. I hope this is a step in the right direction, there seem to be several issues here.

comment:4 Changed 11 years ago by Kyle Shannon

Please let me know how if this makes sense. I think it is a kind of gray area in the driver as, the sphere for the projection is unknown.

comment:5 Changed 10 years ago by etourigny

Keywords: netcdf projection added
Resolution: invalid
Status: assignedclosed

The fundamental problem is that you do not have sufficient information about the existing projection (i.e. the sphere radius). You could encode these using standard CF grid mapping attributes, which are supported in gdal-1.8 and 1.9. Or you could add them in your gdalwarp command using -s_srs '<proj4_string>'.

Also, your target SRS definition (EPSG:3726, valid from -72.0000 to -66.0 west longtitude) does not support the span of your dataset, you should use zone 29.

The following command works (assuming the WGS84 datum):

gdalwarp -s_srs EPSG:4326 -t_srs '+proj=utm +zone=29 +datum=nad83 '

Closing this bug as invalid, because there is nothing wrong with GDAL or the netcdf driver.

Note: See TracTickets for help on using tickets.