[gdal-dev] problem with a geotiff

[Even Rouault]

On mardi 25 septembre 2018 18:59:54 CEST mauge.philippe at free.fr wrote:
> Please, I'd like to submit you an isssue about a geotiff file.
> Go to the web site
> http://www.digitalglobe.com/resources/product-samples/tripoli-lybya Choose
> to download the second file: System-Ready (1B), 50 cm
> Unzip it and try to load the file
> 16MAR08101212-P1BS-055675519010_01_P001.TIF in GDAL I can't load it in GDAL
> because no projection is returned by the fonction GetProjectionRef() It can
> be loaded only if you create manually a .prj and a .tfw file. The file is
> OK because it is correctly loaded in Global Mapper.
> Usually I can load geotiff files in GDAL, but this file won't.
> Can you try to load it and may be tell me what's wrong with it ?

gdalinfo on the file shows:

Size is 26402, 27612
Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
GCP[  0]: Id=1, Info=
          (0.5,0.5) -> (13.156950428616,32.925237256345,0)
GCP[  1]: Id=2, Info=
          (26401.5,0.5) -> (13.298264195371,32.9244661660205,0)
GCP[  2]: Id=3, Info=
          (26401.5,27611.5) -> (13.2982326391638,32.799296089494,0)
GCP[  3]: Id=4, Info=
          (0.5,27611.5) -> (13.157053933887,32.80004957666,0)
Metadata:
  AREA_OR_POINT=Area
  METADATATYPE=DG
  TIFFTAG_COPYRIGHT=(C) COPYRIGHT 2016 DigitalGlobe, Inc., Longmont CO USA 
80503
  TIFFTAG_DATETIME=2016:09:02 06:26:33
  TIFFTAG_IMAGEDESCRIPTION={
  bandList = 
  [
    1;
  ]
}
  TIFFTAG_MAXSAMPLEVALUE=2047
  TIFFTAG_MINSAMPLEVALUE=1
Image Structure Metadata:
  INTERLEAVE=BAND
RPC Metadata:
  HEIGHT_OFF=66
  HEIGHT_SCALE=501
  LAT_OFF=32.8622
  LAT_SCALE=0.0635
  LINE_DEN_COEFF=+1.000000E+00 -4.312485E-05 +9.397242E-04 +3.402903E-05 
+4.831666E-06 +2.711997E-07 -4.030012E-06 -4.962579E-05 +1.425424E-04 
-6.533418E-05 -6.474837E-08 +1.547331E-08 +4.154660E-06 +0.000000E+00 
+4.571326E-07 -2.748659E-04 -5.347683E-08 +0.000000E+00 +5.476937E-06 
+0.000000E+00
  LINE_NUM_COEFF=+1.218335E-03 -6.112947E-03 -1.014484E+00 +7.225448E-03 
-4.622733E-05 +5.629116E-07 +4.233303E-05 -3.719319E-04 -3.608017E-04 
-8.698413E-07 -3.881640E-07 +2.884197E-07 +6.565544E-06 +3.993115E-07 
+4.949166E-05 +2.241984E-04 +6.628490E-05 -3.381833E-07 -2.958911E-06 
-4.712746E-07
  LINE_OFF=13805
  LINE_SCALE=13806
  LONG_OFF=13.2276
  LONG_SCALE=0.0710
  SAMP_DEN_COEFF=+1.000000E+00 -1.207666E-04 +2.303233E-04 -4.443465E-04 
-1.449359E-05 +5.345952E-07 -5.601849E-07 +1.428815E-06 +3.414761E-06 
-2.187619E-06 -2.923988E-08 -1.931437E-08 -2.493269E-08 +0.000000E+00 
+3.951025E-08 -8.725792E-08 +0.000000E+00 +0.000000E+00 +1.117937E-08 
+0.000000E+00
  SAMP_NUM_COEFF=+9.636456E-04 +1.005342E+00 +3.651704E-04 +4.863650E-03 
-2.278309E-04 +4.374349E-04 -2.253067E-04 -8.455434E-04 -5.161798E-04 
+3.445108E-06 +9.375961E-08 -3.540851E-06 -5.042296E-06 -1.878383E-06 
+2.016642E-05 +9.306878E-05 -5.979850E-08 -8.470211E-07 -1.672072E-06 
+0.000000E+00
  SAMP_OFF=13200
  SAMP_SCALE=13201

so the file isn't orthorectified, which explains that GetProjectionRef() 
returns an empty string.

But it contains ground control points (GCPs). You can get the GCP projection 
with GetGCPProjection(). You can compute an approximate geotransform from the 
GCPs with GDALGCPsToGeoTransform(), fed with GDALGetGCPCount() and 
GDALGetGCPs().

Another option is to use the RPC with gdalwarp -rpc (possibly with a DEM with 
-to RPC_DEM=) to orthorectify the file.

Even

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