Opened 2 months ago

Closed 2 months ago

Last modified 2 months ago

#6995 closed defect (fixed)

projection

Reported by: michaeladamkatz Owned by: warmerdam
Priority: normal Milestone: 2.2.2
Component: GDAL_Raster Version: 2.1.3
Severity: normal Keywords: GDALGetProjectionRef
Cc:

Description

In GDAL 2.1.0 and 2.1.3, calling (c++) GDALGetProjectionRef( hDS ) returns empty string for a certain BSB/.kap file (attached). This same code/call worked in GDAL 1.11.

Attachments (1)

18424_2.KAP (496.6 KB) - added by michaeladamkatz 2 months ago.

Download all attachments as: .zip

Change History (16)

Changed 2 months ago by michaeladamkatz

Attachment: 18424_2.KAP added

comment:1 Changed 2 months ago by Even Rouault

In 39788:

GDALGCPsToGeoTransform(): add GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK=YES and GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD=threshold_in_pixel configuration option (refs #6995)

comment:2 Changed 2 months ago by Even Rouault

In 39789:

GDALGCPsToGeoTransform(): add GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK=YES and GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD=threshold_in_pixel configuration option (refs #6995)

comment:3 Changed 2 months ago by Even Rouault

Milestone: 2.2.2
Resolution: fixed
Status: newclosed

The reason is that between GDAL 1.11 and GDAL 2.0 there was a change in the way the GDALGCPsToGeoTransform() function checks the error of compute geotransform from GCP. With the new method, on this dataset, the error is 0.40 pixel whereas the default threshold is 0.25.

$ gdalinfo ~/gdal/data/bsb/18424_2.KAP --debug on
[...]
GDAL: dfErrorX/dfPixelSize = 0.38, dfErrorY/dfPixelSize = 0.40
[...

If you specify the GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD to any value > 0.40 (0.5 for example),

$ gdalinfo ~/gdal/data/bsb/18424_2.KAP --config GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD 0.5
Driver: BSB/Maptech BSB Nautical Charts
Files: /home/even/gdal/data/bsb/18424_2.KAP
Size is 2596, 4652
Coordinate System is:
PROJCS["Global Mercator",
    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"]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-122],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
GeoTransform =
  -58758.74932961759, 1.922511414866929, -6.747722166742181e-05
  6204025.719284005, -2.059470857546899e-05, -1.922518113311302

You could also gdalwarp the .kap into a GeoTIFF to get an orthorectified image.

comment:4 Changed 2 months ago by Even Rouault

Question asked by Michael by email

(1) I'm not only trying to call GDALGetProjectionRef() for this particular
file, but for any file that my program is asked to open. Is there some
number I should set GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD to that will
work for most files? I don't understand what this value does, but I
generally want GDALGetProjectionRef() to return a valid reference string
whenever it can, even at the cost of slight inaccuracy. So I don't know if
I should set it to 1.0, 2.0, or what.

(2) How do I set GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD from the c++
API?

comment:5 Changed 2 months ago by Even Rouault

Michael: I think you can comment on closed tickets, can't you ? (or perhaps it is just I've extended right)

(1) A .KAP file has a set of GCPs that map from pixel coordinates to georeferenced coordinates. From this set, GDALGCPsToGeoTransform() computes a geotransform matrix (affine transformation) used least square adjustments, and then use it to check that the input GCPs are transformed not too far from the georeferenced coordinates they store. This error is measured in pixels. And this is what GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD is about. If you don't care about the inaccuracy, you can just define GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK=YES to disable this post check.

(2) CPLSetConfigOption('GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK', 'YES');

comment:6 Changed 2 months ago by michaeladamkatz

Even, thanks for the explanation.

comment:7 Changed 2 months ago by michaeladamkatz

Unfortunately, adding the call

CPLSetConfigOption( "GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK", "YES" );

either right before or right after the call to

GDALAllRegister();

seems to make no difference.

Is there some other place I should make this call? Some slightly different strings?

comment:8 Changed 2 months ago by michaeladamkatz

Resolution: fixed
Status: closedreopened

comment:9 Changed 2 months ago by Even Rouault

Did you recompile trunk or branches/2.2 from source with the above fixes ?

Last edited 2 months ago by Even Rouault (previous) (diff)

comment:10 Changed 2 months ago by michaeladamkatz

I thought this was a runtime call that would change the behavior dynamically. I am using the precompiled binaries from gisinternals.

comment:11 Changed 2 months ago by Even Rouault

Resolution: fixed
Status: reopenedclosed

This is a runtime call, but you need to grab updated binaries. For example those at http://gisinternals.com/development.php

comment:12 Changed 2 months ago by michaeladamkatz

Well, I updated using the development link you gave. Now I am able to load these .kap files, but now my .sid files are not opening. :(

comment:13 Changed 2 months ago by michaeladamkatz

I have reverted to the stable release of 2.0.3. If there is any other way I can disable the geotransform point checks without going to the development build, please let me know. (Or if you know a quick fix to the .sid problem with the development build, please let me know. But even with that I'd be hesitant to use the development build since other stuff might be broken?)

comment:14 Changed 2 months ago by michaeladamkatz

Bump.

comment:15 Changed 2 months ago by Even Rouault

You need to install the -mrsid.msi packages

Note: See TracTickets for help on using tickets.