[gdal-dev] gdalwarp and geogrids
Even Rouault
even.rouault at spatialys.com
Thu Jan 29 01:12:03 PST 2015
Le jeudi 29 janvier 2015 00:33:48, Even Rouault a écrit :
> 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
Waking up this morning, I realized that the above is necessary, but not sufficient
to help you for your ultimate goal.
gdalwarp will run, but not produce the expected result since it only uses 2D components
of the coordinates. So it will give the same result with and without +geoidgrids
AFAIK, there's no one-step utility to do what you want in GDAL, but you could likely
use the following multiple step process :
1) Convert your grid into a point cloud with gdal2xyz
gdal2xyz.py in.tif > points.txt
2) Use gdaltransform to reproject the point cloud from utm4/navd88 to wgs84
gdaltransform -s_srs "+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs +geoidgrids=./usa_geoid2009/g2009alaska_east_dateline.gtx" \
-t_srs "+proj=longlat +datum=WGS84 +units=m +no_def" < points.txt > points2.txt
3) Edit points2.txt to replace spaces by commas and generate a CSV file
sed "s/ /,/g" < points2.txt > points2.csv
And add a "x,y,z" first line
4) Make a VRT from the CSV to build geometry from the columns. See "Example: ODBC Point Layer" of http://gdal.org/drv_vrt.html as a starting point
5) Optionally convert with ogr2ogr the VRT into a shapefile and add a spatial index for better performance of next step
6) Run gdal_grid to build the output DEM from the VRT/SHP
Even
--
Spatialys - Geospatial professional services
http://www.spatialys.com
More information about the gdal-dev
mailing list