Opened 7 years ago

Closed 7 years ago

Last modified 7 years 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 7 years ago.

Download all attachments as: .zip

Change History (16)

by michaeladamkatz, 7 years ago

Attachment: 18424_2.KAP added

comment:1 by Even Rouault, 7 years ago

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 by Even Rouault, 7 years ago

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 by Even Rouault, 7 years ago

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 by Even Rouault, 7 years ago

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 by Even Rouault, 7 years ago

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 by michaeladamkatz, 7 years ago

Even, thanks for the explanation.

comment:7 by michaeladamkatz, 7 years ago

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 by michaeladamkatz, 7 years ago

Resolution: fixed
Status: closedreopened

comment:9 by Even Rouault, 7 years ago

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

Version 0, edited 7 years ago by Even Rouault (next)

comment:10 by michaeladamkatz, 7 years ago

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

comment:11 by Even Rouault, 7 years ago

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 by michaeladamkatz, 7 years ago

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 by michaeladamkatz, 7 years ago

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 by michaeladamkatz, 7 years ago

Bump.

comment:15 by Even Rouault, 7 years ago

You need to install the -mrsid.msi packages

Note: See TracTickets for help on using tickets.