[gdal-dev] gdalwarp error 6 issues

[Even Rouault]

Le jeudi 01 septembre 2016 19:41:18, Zoltan Szecsei a écrit :
> Hi Even,
> Long time no chat - Hope all is well with you. :-)
> 
> When compiling GDAL, I just used all defaults, so it was just a
> "straight" ./configure.
> 
> The errant images are ENU on WGS84 Lo17 and Lo21, and I need to
> reproject to same but Lo19
> (They are on the fringes of Lo19)
> 
> I tried gdal_edit.py to force a better representation of where those
> images are, but no luck.
> I used variations of the s-srs line below (ie: replacing HBK94 with
> WGS84) and also 1 or 2 EPSG codes - but I was guessing the EPSG ones.
> 
> Is there not a way to stamp out the existing proj info, and replace with
> the desired source projection stanza using gdal_edit.py?

"gdal_edit.py -a_srs srs_def foo.tif" where srs_def is either a proj.4 string 
or EPSG:xxxx code. I've just tried it on a GeoTIFF with a LOCAL_CS originally 
and this worked.

> 
> I've got a batch run in at the moment so don't want to fiddle until that
> has finished.
> I will also recompile with --with-static-proj4.
> 
> Thanks for the pointers.
> Regards,
> Zoltan
> 
> On 2016-09-01 18:40, Even Rouault wrote:
> > Zoltan,
> > 
> >> Hi,
> >> I've been struggling to get some 3200 orthos mosaicked and tiled.
> >> They are very old and have really bad projection definitions, pixel
> >> sizes and metadata funnies.
> >> I'm down to the last 261, and have hit a brick wall.
> >> Given that gdalwarp works on this image when I leave out s_srs and t_srs
> >> (see below), I'm not sure if this is a bug related hiccough or a proj
> >> issue.
> >> 
> >> I have compiled the latest "GDAL 2.1.1, released 2016/07/07" myself, on
> >> Ubuntu 16.04.1 Desktop:  Linux gl05 4.4.0-36-generic #55-Ubuntu SMP Thu
> >> Aug 11 18:01:55 UTC 2016 x86_64 x86_64 x86_64 GNU/Linux
> >> 
> >> I get the following error from gdalwarp:
> >>      ERROR 6: Unable to load PROJ.4 library (libproj.so), creation of
> >>      OGRCoordinateTransformation failed.
> > 
> > Hum, if you get this message, then reprojection should never succeed. How
> > did you compile GDAL with ./configure --with-static-proj4 or without it
> > ? If without it, you must be sure to have the development package of
> > proj installed so that proj.so links to proj.so.1 or proj.so.9 or
> > whatever the soname of your proj is.
> > 
> >> I have tried various different gdalwarp parameters - even leaving
> >> gdalwarp to sense in input projection and me only specifying the output
> >> projection.
> >> If I leave out both s_srs and t_srs then gdalwarp does work but
> >> obliviously does not reproject to Lo19.
> >> I have even tried feeding the output of this redundant gdalwarp into a
> >> gdalwarp with t_srs, but no success.
> >> 
> >> Specifying only t_srs (as above) on many other orthos on Lo17 & Lo21
> >> does work so it is likely that the way this projection is embedded in
> >> this ortho is the problem.
> > 
> > Yes, the way ER Mapper encodes the projection is non-standard and thus
> > GDAL reports a LOCAL_CS that cannot be used to automatically get a
> > proj.4 string. GDAL has some handling for some of those non-standard
> > encodings, but not for your ones, so overriding it with -s_srs is
> > needed.
> > 
> >> I have tried stripping it out with gdal_edit.py but haven't gotten that
> >> right yet.
> >> 
> >> Anyone with pointers or ideas on this?
> > 
> > +datum=Hartebeesthoek94 is not recognized by proj.4.
> > 
> > Why not using EPSG:2049 instead or the proj.4 expension ? I also see that
> > gdalsrsinfo EPSG:2049 expands to
> > 
> > PROJ.4 : '+proj=tmerc +lat_0=0 +lon_0=21 +k=1 +x_0=0 +y_0=0 +axis=wsu
> > +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '
> > 
> > So note the +axis=wsu parameter for the south-orientated variant
> > 
> > Or perhaps this doesn't make sense in your context of producing a
> > GeoTIFF, so in that case you might just want to use the above proj.4
> > string without +axis=wsu. Note that for proj the Hartebeesthoek94 datum
> > is equivalent to WGS84.
> > 
> > Even
> > 
> >> My last effort before email this list is below, as is the gdalinfo
> >> output.
> >> 
> >> Thanks in advance,
> >> Zoltan
> >> 
> >> gdalwarp -s_srs '+proj=tmerc +lat_0=0 +lon_0=21  +k=1 +x_0=0 +y_0=0
> >> +datum=Hartebeesthoek94 +units=m +no_defs'  \
> >> 
> >>            -t_srs '+proj=tmerc +lat_0=0 +lon_0=19  +k=1 +x_0=0 +y_0=0
> >> 
> >> +datum=WGS84 +units=m +no_defs' \
> >> 
> >>            -r cubicspline -multi -overwrite -of GTiff -co TFW=YES -co
> >> 
> >> COMPRESS=JPEG -co JPEG_QUALITY=75 -tr 0.75 -0.75 \
> >> 
> >>            /mnt/sdb1/used_images/dwa_2000/3220/3220CC01.tif \
> >> 
> >> /mnt/sda1/gisdata/raster/used_2016/dwa_2000/3220/3220CC01.tif
> >> #
> >> #
> >> gdalinfo   /mnt/sdb1/used_images/dwa_2000/3220/3220CC01.tif
> >> #
> >> Driver: GTiff/GeoTIFF
> >> Files: /mnt/sdb1/used_images/dwa_2000/3220/3220CC01.tif
> >> Size is 6729, 7720
> >> Coordinate System is:
> >> LOCAL_CS[" Geocoding information not available Projection Name = WGS
> >> 1984 Lo21 Units = meters GeoTIFF Units = meters",
> >> 
> >>       UNIT["metre",1,
> >>       
> >>           AUTHORITY["EPSG","9001"]]]
> >> 
> >> Origin = (-93872.189308887900552,-3625321.216973003931344)
> >> Pixel Size = (0.750000000000000,-0.750000000000000)
> >> 
> >> Metadata:
> >>     AREA_OR_POINT=Area
> >>     TIFFTAG_DOCUMENTNAME=hetotiff
> >>     TIFFTAG_IMAGEDESCRIPTION=ER Mapper 6.4 hetotiff
> >>     TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
> >>     TIFFTAG_SOFTWARE=IMAGINE TIFF Support
> >> 
> >> Copyright 1991 - 1999 by ERDAS, Inc. All Rights Reserved
> >> @(#)$RCSfile: etif.c $ $Revision: 1.9.1.3 $ $Date: 2002/07/29
> >> 15:39:06EDT $
> >> 
> >>     TIFFTAG_XRESOLUTION=1
> >>     TIFFTAG_YRESOLUTION=1
> >> 
> >> Image Structure Metadata:
> >>     INTERLEAVE=PIXEL
> >> 
> >> Corner Coordinates:
> >> Upper Left  (  -93872.189,-3625321.217)
> >> Lower Left  (  -93872.189,-3631111.217)
> >> Upper Right (  -88825.439,-3625321.217)
> >> Lower Right (  -88825.439,-3631111.217)
> >> Center      (  -91348.814,-3628216.217)
> >> Band 1 Block=6729x1 Type=Byte, ColorInterp=Red
> >> Band 2 Block=6729x1 Type=Byte, ColorInterp=Green
> >> Band 3 Block=6729x1 Type=Byte, ColorInterp=Blue

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