Error while reprojecting with gdalwarp an ECMWF netCDF

[toyzerocha]

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

[toyzerocha]

NetCDF original file

[toyzerocha]

Tif file after gdalwarp

[kyle]

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.

[kyle]

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.

[etourigny]

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.