[gdal-dev] gdalwarp and geogrids

[Even Rouault]

Le mercredi 28 janvier 2015 21:53:12, Wood, Alexander a écrit :
> Hello,
> 
> Hello devlist,
> I have small patch of alaska elevation data  in utm4/navd88 and I would
> like to convert this into wgs84 ellipsoidal.  I know this data uses
> navd88,but gdalinfo does not report this vertical datum in the WKT.  To
> reproject this data into wgs84 ellipsoidal elevation, I'm assuming it
> should be possible with the following gdalwarp command:
> 
> gdalwarp -s_srs "+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs
> +geoidgrids=g2009alaska.gtx" -t_srs "+proj=longlat +datum=WGS84 +units=m
> +no_def" alaska_sample_utm4_navd88.tif alaska_sample_wgs84.tif
> 
> I'm seeing the following errors reported from this command:
> ERROR 1: point not within available datum shift grids
> ERROR 1: GDALWarperOperation::ComputeSourceWindow() failed because the
> pfnTransformer failed. Warning 1: Unable to compute source region for
> output window 0,0,351,106, skipping.
> 
> I'm using the geoidgrids from
> http://trac.osgeo.org/proj/wiki/VerticalDatums, moved into the
> gdal/projlib directory
> 
> The source file is in utm4:
> PROJ.4 string is:
> '+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs '
> Origin = (412011.848695503660000,7085064.815277445100000)
> 
> For reference, this is about 100msurrounding the following lat/long:
> Center      (-160.7909106,  63.8815138) (160d47'27.28"W, 63d52'53.45"N)
> 
> A simple gdalinfo on g2009alaska.gtx reports the following:
> PROJ.4 string is:
> '+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs '
> Origin = (171.991666666666670,72.008333333333340)
> Pixel Size = (0.016666666666667,-0.016666666666667)
> Corner Coordinates:
> Upper Left  ( 171.9916667,  72.0083333) (171d59'30.00"E, 72d 0'30.00"N)
> Lower Left  ( 171.9916667,  48.9916667) (171d59'30.00"E, 48d59'30.00"N)
> Upper Right (     234.008,      72.008) (234d 0'30.00"E, 72d 0'30.00"N)
> Lower Right (     234.008,      48.992) (234d 0'30.00"E, 48d59'30.00"N)
> Center      (     203.000,      60.500) (203d 0' 0.00"E, 60d30' 0.00"N)
> 
> Does anyone see anything wrong with my gdalwarp/proj4 usage?  Are there
> better ways to handle these sorts of transformations?

Alexander,

I can also reproduce. Running with PROJ_DEBUG=ON gives an interesting hint :

$ echo "-160.7909106 63.8815138" | PROJ_DEBUG=ON gdaltransform  -s_srs 
"+proj=longlat +datum=WGS84 +units=m +no_defs 
+geoidgrids=./usa_geoid2009/g2009alaska.gtx" -t_srs "+proj=longlat 
+datum=WGS84 +units=m +no_def" --debug on
OGRCT: PROJ >= 4.8.0 features enabled
pj_open_lib(proj_def.dat): call fopen(/home/even/install-
proj-4.8.0/share/proj/proj_def.dat) - succeeded

OGRCT: Source: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 
+geoidgrids=./usa_geoid2009/g2009alaska.gtx +no_defs 
OGRCT: Target: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
OGRCT: Source: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
OGRCT: Target: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 
+geoidgrids=./usa_geoid2009/g2009alaska.gtx +no_defs 
pj_open_lib(./usa_geoid2009/g2009alaska.gtx): call 
fopen(./usa_geoid2009/g2009alaska.gtx) - succeeded

This GTX spans the dateline!  This will cause problems.
GTX 3721x1381: LL=(172,49) UR=(234,72)
pj_apply_vgridshift(): failed to find a grid shift table for
                       location (-160.7909106dW,63.8815138dN)
   tried: ./usa_geoid2009/g2009alaska.gtx
ERROR 1: point not within available datum shift grids
transformation failed.

And indeed :
$ gdalinfo ./usa_geoid2009/g2009alaska.gtx
Driver: GTX/NOAA Vertical Datum .GTX
Files: ./usa_geoid2009/g2009alaska.gtx
Size is 3721, 1381
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 = (171.991666666666674,72.008333333333340)
Pixel Size = (0.016666666666667,-0.016666666666667)
Corner Coordinates:
Upper Left  ( 171.9916667,  72.0083333) (171d59'30.00"E, 72d 0'30.00"N)
Lower Left  ( 171.9916667,  48.9916667) (171d59'30.00"E, 48d59'30.00"N)
Upper Right (     234.008,      72.008) (234d 0'30.00"E, 72d 0'30.00"N)
Lower Right (     234.008,      48.992) (234d 0'30.00"E, 48d59'30.00"N)
Center      (     203.000,      60.500) (203d 0' 0.00"E, 60d30' 0.00"N)
Band 1 Block=3721x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-88.8888015747070312

Current proj.4 code can deal with grids whose longitudes are all < 180 or all 
> 180, but not crossing. This could probably be improved, but would require 
further investigation.

I've successfully tried the following workaround by cutting the grid into 2 
parts around the dateline :

$ gdal_translate ./usa_geoid2009/g2009alaska.gtx \
			 ./usa_geoid2009/g2009alaska_west_dateline.gtx -of gtx \
			 -projwin 171.9916667 72.0083333 180 48.992


$ gdal_translate ./usa_geoid2009/g2009alaska.gtx \
			 ./usa_geoid2009/g2009alaska_east_dateline.gtx -of gtx \
			 -projwin 180 72.0083333 234.008 48.992


And then :
$ echo "-160.7909106 63.8815138" | gdaltransform  -s_srs "+proj=longlat 
+datum=WGS84 +units=m +no_defs 
+geoidgrids=./usa_geoid2009/g2009alaska_east_dateline.gtx" -t_srs 
"+proj=longlat +datum=WGS84 +units=m +no_def"
-160.7909106 63.8815138 7.7247364547829

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com