Skip to content

Commit

Permalink
GTiff: only report scale/offset deduced from ModelTiepointTag and Mod…
Browse files Browse the repository at this point in the history 
…elPixelScaleTag if the SRS has a vertical component (and thus currently if GTIFF_REPORT_COMPD_CS is set) (fixes https://issues.qgis.org/issues/20493)
  • Loading branch information
@rouault
rouault committed Jan 30, 2019
1 parent b412b0e commit 0a3d241
Show file tree
Hide file tree
Showing 4 changed files with 74 additions and 57 deletions.
Binary file added autotest/gcore/data/ModelTiepointTag_z_non_zero_but_ModelPixelScaleTag_z_zero.tif
View file
Binary file not shown.
8 changes: 6 additions & 2 deletions autotest/gcore/tiff_read.py
View file
Expand Up @@ -3119,7 +3119,11 @@ def test_tiff_read_overview_of_external_mask():
assert flags1 == gdal.GMF_PER_DATASET

###############################################################################
# Test reading GeoTIFF file ModelTiepointTag(z) != 0 and ModelPixelScaleTag(z) = 0
# Test https://issues.qgis.org/issues/20493

def test_tiff_read_ModelTiepointTag_z_non_zero_but_ModelPixelScaleTag_z_zero():



ds = gdal.Open('data/ModelTiepointTag_z_non_zero_but_ModelPixelScaleTag_z_zero.tif')
assert ds.GetRasterBand(1).GetScale() == 1
assert ds.GetRasterBand(1).GetOffset() == 0
87 changes: 45 additions & 42 deletions autotest/gcore/tiff_write.py
View file
Expand Up @@ -3557,33 +3557,31 @@ def test_tiff_write_97():

def test_tiff_write_98():

gdal.SetConfigOption('GTIFF_POINT_GEO_IGNORE', 'FALSE')

src_ds = gdal.Open('data/geomatrix.tif')
with gdaltest.config_option('GTIFF_POINT_GEO_IGNORE', 'FALSE'):
src_ds = gdal.Open('data/geomatrix.tif')

gdal.SetConfigOption('GTIFF_POINT_GEO_IGNORE', 'TRUE')
with gdaltest.config_option('GTIFF_POINT_GEO_IGNORE', 'TRUE'):
new_ds = gdaltest.tiff_drv.CreateCopy('tmp/test_98.tif', src_ds)

new_ds = gdaltest.tiff_drv.CreateCopy('tmp/test_98.tif', src_ds)

gt = new_ds.GetGeoTransform()
md = new_ds.GetMetadataItem('AREA_OR_POINT')
new_ds = None
src_ds = None
gt = new_ds.GetGeoTransform()
md = new_ds.GetMetadataItem('AREA_OR_POINT')
new_ds = None
src_ds = None

gt_expected = (1841001.75, 1.5, -5.0, 1144003.25, -5.0, -1.5)
gt_expected = (1841001.75, 1.5, -5.0, 1144003.25, -5.0, -1.5)

assert gt == gt_expected, 'did not get expected geotransform'

assert md == 'Point', 'did not get expected AREA_OR_POINT value'

gdal.SetConfigOption('GTIFF_POINT_GEO_IGNORE', 'FALSE')
with gdaltest.config_option('GTIFF_POINT_GEO_IGNORE', 'FALSE'):

new_ds = gdal.Open('tmp/test_98.tif')
new_ds = gdal.Open('tmp/test_98.tif')

gt = new_ds.GetGeoTransform()
md = new_ds.GetMetadataItem('AREA_OR_POINT')
new_ds = None
src_ds = None
gt = new_ds.GetGeoTransform()
md = new_ds.GetMetadataItem('AREA_OR_POINT')
new_ds = None
src_ds = None

gt_expected = (1841003.5, 1.5, -5.0, 1144006.5, -5.0, -1.5)

Expand Down Expand Up @@ -3727,14 +3725,14 @@ def test_tiff_write_102():
ds.SetProjection(wkt)
ds = None

gdal.SetConfigOption('GTIFF_REPORT_COMPD_CS', 'YES')
ds = gdal.Open('/vsimem/tiff_write_102.tif')
wkt1 = ds.GetProjectionRef()
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('/vsimem/tiff_write_102.tif')
wkt1 = ds.GetProjectionRef()
ds = None

gdal.SetConfigOption('GTIFF_REPORT_COMPD_CS', 'NO')
ds = gdal.Open('/vsimem/tiff_write_102.tif')
wkt2 = ds.GetProjectionRef()
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'NO'):
ds = gdal.Open('/vsimem/tiff_write_102.tif')
wkt2 = ds.GetProjectionRef()
ds = None

gdaltest.tiff_drv.Delete('/vsimem/tiff_write_102.tif')
Expand Down Expand Up @@ -6362,20 +6360,21 @@ def test_tiff_write_165():

def test_tiff_write_166():

ds = gdal.Open('data/tiff_vertcs_scale_offset.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0

assert ds.GetRasterBand(1).GetOffset() == 10.0
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('data/tiff_vertcs_scale_offset.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
assert ds.GetRasterBand(1).GetOffset() == 10.0

# Scale + offset through CreateCopy()
gdal.Translate('/vsimem/tiff_write_166.tif', 'data/byte.tif',
options='-a_srs EPSG:26711+5773 -a_scale 2.0 -a_offset 10 -co PROFILE=GEOTIFF')
assert gdal.VSIStatL('/vsimem/tiff_write_166.tif.aux.xml') is None

ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
assert ds.GetRasterBand(1).GetOffset() == 10.0

assert ds.GetRasterBand(1).GetOffset() == 10.0
ds = None
gdal.Unlink('/vsimem/tiff_write_166.tif')

Expand All @@ -6384,10 +6383,11 @@ def test_tiff_write_166():
options='-a_srs EPSG:26711+5773 -a_offset 10 -co PROFILE=GEOTIFF')
assert gdal.VSIStatL('/vsimem/tiff_write_166.tif.aux.xml') is None

ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 1.0
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 1.0
assert ds.GetRasterBand(1).GetOffset() == 10.0

assert ds.GetRasterBand(1).GetOffset() == 10.0
ds = None
gdal.Unlink('/vsimem/tiff_write_166.tif')

Expand All @@ -6402,9 +6402,10 @@ def test_tiff_write_166():
ds = None
assert gdal.VSIStatL('/vsimem/tiff_write_166.tif.aux.xml') is None

ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
assert ds.GetRasterBand(1).GetOffset() == 10.0
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
assert ds.GetRasterBand(1).GetOffset() == 10.0
ds = None
gdal.Unlink('/vsimem/tiff_write_166.tif')

Expand All @@ -6418,9 +6419,10 @@ def test_tiff_write_166():
ds = None
assert gdal.VSIStatL('/vsimem/tiff_write_166.tif.aux.xml') is None

ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
assert ds.GetRasterBand(1).GetOffset() == 0.0
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 2.0
assert ds.GetRasterBand(1).GetOffset() == 0.0
ds = None
gdal.Unlink('/vsimem/tiff_write_166.tif')

Expand All @@ -6434,9 +6436,10 @@ def test_tiff_write_166():
ds = None
assert gdal.VSIStatL('/vsimem/tiff_write_166.tif.aux.xml') is None

ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 1.0
assert ds.GetRasterBand(1).GetOffset() == 10.0
with gdaltest.config_option('GTIFF_REPORT_COMPD_CS', 'YES'):
ds = gdal.Open('/vsimem/tiff_write_166.tif')
assert ds.GetRasterBand(1).GetScale() == 1.0
assert ds.GetRasterBand(1).GetOffset() == 10.0
ds = None
gdal.Unlink('/vsimem/tiff_write_166.tif')

Expand Down
36 changes: 23 additions & 13 deletions gdal/frmts/gtiff/geotiff.cpp
View file
Expand Up @@ -14615,19 +14615,23 @@ void GTiffDataset::LoadGeoreferencingAndPamIfNeeded()
padfTiePoints[2] != 0.0 ||
padfTiePoints[5] != 0.0) )
{
/* modelTiePointTag = (pixel, line, z0, X, Y, Z0) */
/* thus Z(some_point) = (z(some_point) - z0) * scaleZ + Z0 */
/* equivalently written as */
/* Z(some_point) = z(some_point) * scaleZ + offsetZ with */
/* offsetZ = - z0 * scaleZ + Z0 */
double dfScale = padfScale[2];
double dfOffset =
-padfTiePoints[2] * dfScale + padfTiePoints[5];
GTiffRasterBand* poBand =
cpl::down_cast<GTiffRasterBand*>(GetRasterBand(1));
poBand->bHaveOffsetScale = true;
poBand->dfScale = dfScale;
poBand->dfOffset = dfOffset;
LookForProjection();
if( pszProjection && HasVerticalCS(pszProjection) )
{
/* modelTiePointTag = (pixel, line, z0, X, Y, Z0) */
/* thus Z(some_point) = (z(some_point) - z0) * scaleZ + Z0 */
/* equivalently written as */
/* Z(some_point) = z(some_point) * scaleZ + offsetZ with */
/* offsetZ = - z0 * scaleZ + Z0 */
double dfScale = padfScale[2];
double dfOffset =
-padfTiePoints[2] * dfScale + padfTiePoints[5];
GTiffRasterBand* poBand =
cpl::down_cast<GTiffRasterBand*>(GetRasterBand(1));
poBand->bHaveOffsetScale = true;
poBand->dfScale = dfScale;
poBand->dfOffset = dfOffset;
}
}
}
}
Expand Down Expand Up @@ -17459,7 +17463,13 @@ GTiffDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
}
}

CPLString osOldGTIFF_REPORT_COMPD_CSVal(
CPLGetConfigOption("GTIFF_REPORT_COMPD_CS", ""));
CPLSetThreadLocalConfigOption("GTIFF_REPORT_COMPD_CS", "YES");
poDS->CloneInfo( poSrcDS, nCloneInfoFlags );
CPLSetThreadLocalConfigOption("GTIFF_REPORT_COMPD_CS",
osOldGTIFF_REPORT_COMPD_CSVal.empty() ? nullptr :
osOldGTIFF_REPORT_COMPD_CSVal.c_str());

if( !bGeoTIFF && (poDS->GetPamFlags() & GPF_DISABLED) == 0 )
{
Expand Down

0 comments on commit 0a3d241

Please sign in to comment.