Skip to content

Commit

Permalink
Merge pull request OSGeo#7730 from rouault/isg_arg_fix
Browse files Browse the repository at this point in the history 
ISG: relax tolerence check to accept GEOIDEAR16_20160419.isg
  • Loading branch information
@rouault
rouault authored May 23, 2023
2 parents e48a8f8 + e7906a1 commit 65cb355
Show file tree
Hide file tree
Showing 5 changed files with 210 additions and 78 deletions.
16 changes: 16 additions & 0 deletions autotest/gdrivers/data/isg/approx_georeferencing_auto_corrected.isg
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
begin_of_head ================================================
model name : GEOIDEAR16
model type : gravimetric
units : meters
reference : --
lat min = -57.0091
lat max = -20.0083
lon min = -76.0083
lon max = -52.0079
delta lat = 0.0167
delta lon = 0.0167
nrows = 2220
ncols = 1440
nodata = -9999.0000
ISG format = 1.0
end_of_head ==================================================
17 changes: 17 additions & 0 deletions autotest/gdrivers/data/isg/approx_georeferencing_rejected_by_default.isg
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
begin_of_head ================================================
model name : GEOIDEAR16
model type : gravimetric
units : meters
reference : --
lat min = -57.0091
lat max = -20.0083
lon min = -76.0083
lon max = -52.0079
delta lat = 0.0167
# delta lon should be 0.0167
delta lon = 0.0165
nrows = 2220
ncols = 1440
nodata = -9999.0000
ISG format = 1.0
end_of_head ==================================================
17 changes: 17 additions & 0 deletions autotest/gdrivers/data/isg/approx_georeferencing_warning.isg
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
begin_of_head ================================================
model name : GEOIDEAR16
model type : gravimetric
units : meters
reference : --
lat min = -57.0091
lat max = -20.0083
lon min = -76.0083
lon max = -52.0079
delta lat = 0.0167
# delta lon should be 0.0167
delta lon = 0.0166
nrows = 2220
ncols = 1440
nodata = -9999.0000
ISG format = 1.0
end_of_head ==================================================
71 changes: 71 additions & 0 deletions autotest/gdrivers/isg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,11 @@
###############################################################################

import gdaltest
import pytest

from osgeo import gdal

pytestmark = pytest.mark.require_driver("ISG")

###############################################################################
# Perform simple read test.
Expand All @@ -39,3 +44,69 @@ def test_isg_1():
tst = gdaltest.GDALTest("ISG", "isg/test.isg", 1, 159)
expected_gt = [120.0, 0.25, 0.0, 41.0, 0.0, -0.25]
return tst.testOpen(check_gt=expected_gt)


###############################################################################
# Test if we can accept approximate georeferencing


def test_isg_approx_georeferencing_auto_corrected():

gdal.ErrorReset()
# Header of https://www.isgeoid.polimi.it/Geoid/America/Argentina/public/GEOIDEAR16_20160419.isg
ds = gdal.Open("data/isg/approx_georeferencing_auto_corrected.isg")
assert gdal.GetLastErrorMsg() == ""
expected_gt = [-76.0098535, 0.016667, 0.0, -20.0087335, 0.0, -0.016667]
assert ds.GetGeoTransform() == pytest.approx(expected_gt, rel=1e-8)


###############################################################################
# Test if we can accept approximate georeferencing


def test_isg_approx_georeferencing_with_warning():

gdal.ErrorReset()
# Header of https://www.isgeoid.polimi.it/Geoid/America/Argentina/public/GEOIDEAR16_20160419.isg
# with delta_lon modified
with gdaltest.error_handler():
ds = gdal.Open("data/isg/approx_georeferencing_warning.isg")
assert gdal.GetLastErrorMsg() != ""
expected_gt = [
-76.0083,
0.01666694444444445,
0.0,
-20.0083,
0.0,
-0.016667027027027027,
]
assert ds.GetGeoTransform() == pytest.approx(expected_gt, rel=1e-8)


###############################################################################
# Test if we can accept approximate georeferencing


def test_isg_approx_georeferencing_rejected_by_default():

# Header of https://www.isgeoid.polimi.it/Geoid/America/Argentina/public/GEOIDEAR16_20160419.isg
# with delta_lon modified
filename = "data/isg/approx_georeferencing_rejected_by_default.isg"
with pytest.raises(Exception, match="ISG_SKIP_GEOREF_CONSISTENCY_CHECK"):
assert gdal.Open(filename) is None

gdal.ErrorReset()
with gdal.config_option(
"ISG_SKIP_GEOREF_CONSISTENCY_CHECK", "YES"
), gdaltest.error_handler():
ds = gdal.Open(filename)
assert gdal.GetLastErrorMsg() != ""
expected_gt = [
-76.0083,
0.01666694444444445,
0.0,
-20.0083,
0.0,
-0.016667027027027027,
]
assert ds.GetGeoTransform() == pytest.approx(expected_gt, rel=1e-8)
167 changes: 89 additions & 78 deletions frmts/aaigrid/aaigriddataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -823,105 +823,116 @@ int ISGDataset::ParseHeader(const char *pszHeader, const char *)

// Correct rounding errors.

const double dfRoundedDeltaLon =
(osDeltaLon == "0.0167" ||
(dfDeltaLon < 1 &&
fabs(1. / dfDeltaLon - floor(1. / dfDeltaLon + 0.5)) < 0.06))
? 1. / floor(1. / dfDeltaLon + 0.5)
: dfDeltaLon;
if (dfRoundedDeltaLon != dfDeltaLon &&
fabs(fabs(dfLonMin / dfRoundedDeltaLon) -
(floor(fabs(dfLonMin / dfRoundedDeltaLon)) + 0.5)) < 0.02 &&
fabs(fabs(dfLonMax / dfRoundedDeltaLon) -
(floor(fabs(dfLonMax / dfRoundedDeltaLon)) + 0.5)) < 0.02)
const auto TryRoundTo = [](double &dfDelta, double dfRoundedDelta,
double &dfMin, double &dfMax, int nVals,
double dfRelTol)
{
double dfMinTry = dfMin;
double dfMaxTry = dfMax;
double dfDeltaTry = dfDelta;
if (dfRoundedDelta != dfDelta &&
fabs(fabs(dfMin / dfRoundedDelta) -
(floor(fabs(dfMin / dfRoundedDelta)) + 0.5)) < dfRelTol &&
fabs(fabs(dfMax / dfRoundedDelta) -
(floor(fabs(dfMax / dfRoundedDelta)) + 0.5)) < dfRelTol)
{
double dfVal = (floor(fabs(dfLonMin / dfRoundedDeltaLon)) + 0.5) *
dfRoundedDeltaLon;
dfLonMin = (dfLonMin < 0) ? -dfVal : dfVal;
}
{
double dfVal = (floor(fabs(dfLonMax / dfRoundedDeltaLon)) + 0.5) *
dfRoundedDeltaLon;
dfLonMax = (dfLonMax < 0) ? -dfVal : dfVal;
{
double dfVal = (floor(fabs(dfMin / dfRoundedDelta)) + 0.5) *
dfRoundedDelta;
dfMinTry = (dfMin < 0) ? -dfVal : dfVal;
}
{
double dfVal = (floor(fabs(dfMax / dfRoundedDelta)) + 0.5) *
dfRoundedDelta;
dfMaxTry = (dfMax < 0) ? -dfVal : dfVal;
}
dfDeltaTry = dfRoundedDelta;
}
dfDeltaLon = dfRoundedDeltaLon;
}
else if (dfRoundedDeltaLon != dfDeltaLon &&
fabs(fabs(dfLonMin / dfRoundedDeltaLon) -
(floor(fabs(dfLonMin / dfRoundedDeltaLon) + 0.5) + 0.)) <
0.02 &&
fabs(fabs(dfLonMax / dfRoundedDeltaLon) -
(floor(fabs(dfLonMax / dfRoundedDeltaLon) + 0.5) + 0.)) <
0.02)
{
else if (dfRoundedDelta != dfDelta &&
fabs(fabs(dfMin / dfRoundedDelta) -
(floor(fabs(dfMin / dfRoundedDelta) + 0.5) + 0.)) <
dfRelTol &&
fabs(fabs(dfMax / dfRoundedDelta) -
(floor(fabs(dfMax / dfRoundedDelta) + 0.5) + 0.)) <
dfRelTol)
{
double dfVal =
(floor(fabs(dfLonMin / dfRoundedDeltaLon) + 0.5) + 0.) *
dfRoundedDeltaLon;
dfLonMin = (dfLonMin < 0) ? -dfVal : dfVal;
{
double dfVal =
(floor(fabs(dfMin / dfRoundedDelta) + 0.5) + 0.) *
dfRoundedDelta;
dfMinTry = (dfMin < 0) ? -dfVal : dfVal;
}
{
double dfVal =
(floor(fabs(dfMax / dfRoundedDelta) + 0.5) + 0.) *
dfRoundedDelta;
dfMaxTry = (dfMax < 0) ? -dfVal : dfVal;
}
dfDeltaTry = dfRoundedDelta;
}
if (fabs(dfMinTry + dfDeltaTry * nVals - dfMaxTry) <
dfRelTol * dfDeltaTry)
{
double dfVal =
(floor(fabs(dfLonMax / dfRoundedDeltaLon) + 0.5) + 0.) *
dfRoundedDeltaLon;
dfLonMax = (dfLonMax < 0) ? -dfVal : dfVal;
dfMin = dfMinTry;
dfMax = dfMaxTry;
dfDelta = dfDeltaTry;
return true;
}
dfDeltaLon = dfRoundedDeltaLon;
}
return false;
};

const double dfRoundedDeltaLon =
(osDeltaLon == "0.0167" ||
(dfDeltaLon < 1 &&
fabs(1. / dfDeltaLon - floor(1. / dfDeltaLon + 0.5)) < 0.06))
? 1. / floor(1. / dfDeltaLon + 0.5)
: dfDeltaLon;

const double dfRoundedDeltaLat =
(osDeltaLat == "0.0167" ||
(dfDeltaLat < 1 &&
fabs(1. / dfDeltaLat - floor(1. / dfDeltaLat + 0.5)) < 0.06))
? 1. / floor(1. / dfDeltaLat + 0.5)
: dfDeltaLat;
if (dfRoundedDeltaLat != dfDeltaLat &&
fabs(fabs(dfLatMin / dfRoundedDeltaLat) -
(floor(fabs(dfLatMin / dfRoundedDeltaLat)) + 0.5)) < 0.02 &&
fabs(fabs(dfLatMax / dfRoundedDeltaLat) -
(floor(fabs(dfLatMax / dfRoundedDeltaLat)) + 0.5)) < 0.02)

bool bOK = TryRoundTo(dfDeltaLon, dfRoundedDeltaLon, dfLonMin, dfLonMax,
nCols, 1e-2) &&
TryRoundTo(dfDeltaLat, dfRoundedDeltaLat, dfLatMin, dfLatMax,
nRows, 1e-2);
if (!bOK && osDeltaLon == "0.0167" && osDeltaLat == "0.0167")
{
{
double dfVal = (floor(fabs(dfLatMin / dfRoundedDeltaLat)) + 0.5) *
dfRoundedDeltaLat;
dfLatMin = (dfLatMin < 0) ? -dfVal : dfVal;
}
{
double dfVal = (floor(fabs(dfLatMax / dfRoundedDeltaLat)) + 0.5) *
dfRoundedDeltaLat;
dfLatMax = (dfLatMax < 0) ? -dfVal : dfVal;
}
dfDeltaLat = dfRoundedDeltaLat;
// For https://www.isgeoid.polimi.it/Geoid/America/Argentina/public/GEOIDEAR16_20160419.isg
bOK =
TryRoundTo(dfDeltaLon, 0.016667, dfLonMin, dfLonMax, nCols, 1e-1) &&
TryRoundTo(dfDeltaLat, 0.016667, dfLatMin, dfLatMax, nRows, 1e-1);
}
else if (dfRoundedDeltaLat != dfDeltaLat &&
fabs(fabs(dfLatMin / dfRoundedDeltaLat) -
(floor(fabs(dfLatMin / dfRoundedDeltaLat) + 0.5) + 0.)) <
0.02 &&
fabs(fabs(dfLatMax / dfRoundedDeltaLat) -
(floor(fabs(dfLatMax / dfRoundedDeltaLat) + 0.5) + 0.)) <
0.02)
if (!bOK)
{
// 0.005 is what would be needed for the above GEOIDEAR16_20160419.isg
// file without the specific fine tuning done.
if ((fabs((dfLonMax - dfLonMin) / nCols - dfDeltaLon) <
0.005 * dfDeltaLon &&
fabs((dfLatMax - dfLatMin) / nRows - dfDeltaLat) <
0.005 * dfDeltaLat) ||
CPLTestBool(
CPLGetConfigOption("ISG_SKIP_GEOREF_CONSISTENCY_CHECK", "NO")))
{
double dfVal =
(floor(fabs(dfLatMin / dfRoundedDeltaLat) + 0.5) + 0.) *
dfRoundedDeltaLat;
dfLatMin = (dfLatMin < 0) ? -dfVal : dfVal;
CPLError(CE_Warning, CPLE_AppDefined,
"Georeference might be slightly approximate due to "
"rounding of coordinates and resolution in file header.");
dfDeltaLon = (dfLonMax - dfLonMin) / nCols;
dfDeltaLat = (dfLatMax - dfLatMin) / nRows;
}
else
{
double dfVal =
(floor(fabs(dfLatMax / dfRoundedDeltaLat) + 0.5) + 0.) *
dfRoundedDeltaLat;
dfLatMax = (dfLatMax < 0) ? -dfVal : dfVal;
CPLError(CE_Failure, CPLE_AppDefined,
"Inconsistent extent/resolution/raster dimension, or "
"rounding of coordinates and resolution in file header "
"higher than accepted. You may skip this consistency "
"check by setting the ISG_SKIP_GEOREF_CONSISTENCY_CHECK "
"configuration option to YES.");
return false;
}
dfDeltaLat = dfRoundedDeltaLat;
}

if (!(fabs(dfLatMin + dfDeltaLat * nRows - dfLatMax) < 1e-8 &&
fabs(dfLonMin + dfDeltaLon * nCols - dfLonMax) < 1e-8))
{
CPLDebug("ISG", "Inconsistent extent/resolution/raster dimension");
return FALSE;
}
nRasterXSize = nCols;
nRasterYSize = nRows;
Expand Down

0 comments on commit 65cb355

Please sign in to comment.