Skip to content

Commit

Permalink
Merge pull request #3182 from rouault/proj8_merc
Browse files Browse the repository at this point in the history 
gdal2tiles/COG/MBTiles/GeoPackage: adjustments for EPSG:3857 output (due to PROJ 8 changes)
  • Loading branch information
@rouault
rouault committed Nov 15, 2020
2 parents e86cbc4 + fd8615b commit dc38aa6
Show file tree
Hide file tree
Showing 7 changed files with 266 additions and 153 deletions.
4 changes: 3 additions & 1 deletion autotest/gdrivers/mbtiles.py
View file
Expand Up @@ -435,8 +435,10 @@ def test_mbtiles_8():
out_ds = None
src_ds = None

expected_cs = [60844, 7388, 53813]
expected_cs = [580, 8742, 54747]
out_ds = gdal.Open('/vsimem/mbtiles_8.mbtiles')
assert out_ds.RasterXSize == 512
assert out_ds.RasterYSize == 512
got_cs = [out_ds.GetRasterBand(i + 1).Checksum() for i in range(3)]
assert got_cs == expected_cs
got_ct = out_ds.GetRasterBand(1).GetColorTable()
Expand Down
8 changes: 4 additions & 4 deletions autotest/pyscripts/test_gdal2tiles.py
View file
Expand Up @@ -70,7 +70,7 @@ def test_gdal2tiles_py_simple():

_verify_raster_band_checksums(
'tmp/out_gdal2tiles_smallworld/0/0/0.png',
expected_cs = [25314, 28114, 6148, 59026]
expected_cs = [31420, 32522, 16314, 17849]
)

for filename in ['googlemaps.html', 'leaflet.html', 'openlayers.html', 'tilemapresource.xml']:
Expand All @@ -97,7 +97,7 @@ def test_gdal2tiles_py_zoom_option():

_verify_raster_band_checksums(
'tmp/out_gdal2tiles_smallworld/1/0/0.png',
expected_cs = [8130, 10496, 65274, 63715]
expected_cs = [24063, 23632, 14707, 17849]
)

ds = gdal.Open('tmp/out_gdal2tiles_smallworld/doc.kml')
Expand Down Expand Up @@ -164,11 +164,11 @@ def test_gdal2tiles_py_xyz():

_verify_raster_band_checksums(
'tmp/out_gdal2tiles_smallworld_xyz/0/0/0.png',
expected_cs = [30616, 31851, 9392, 63557]
expected_cs = [31747, 33381, 18447, 17849]
)
_verify_raster_band_checksums(
'tmp/out_gdal2tiles_smallworld_xyz/1/0/0.png',
expected_cs = [25095, 27337, 10068, 63699]
expected_cs = [15445, 16942, 13681, 17849]
)

for filename in ['googlemaps.html', 'leaflet.html', 'openlayers.html']:
Expand Down
14 changes: 10 additions & 4 deletions autotest/utilities/test_gdalwarp.py
View file
Expand Up @@ -599,7 +599,11 @@ def test_gdalwarp_28():
# Test warping a full EPSG:4326 extent to EPSG:3785 (#2305)


def test_gdalwarp_29():
def DISABLED_test_gdalwarp_29():

# This test has been disabled since PROJ 8 will reproject a coordinates at
# lat=90 to a finite value, due to 90deg being < PI/2 due to numerical
# accuracy
if test_cli_utilities.get_gdalwarp_path() is None:
pytest.skip()

Expand All @@ -626,14 +630,16 @@ def test_gdalwarp_30():
if test_cli_utilities.get_gdalwarp_path() is None:
pytest.skip()

te = " -te -20037508.343 -16206629.152 20036845.112 16213801.068"

# First run : no parameter
gdaltest.runexternal(test_cli_utilities.get_gdalwarp_path() + " data/w_jpeg.tiff tmp/testgdalwarp30_1.tif -t_srs EPSG:3785 -co COMPRESS=LZW -wm 500000 --config GDAL_CACHEMAX 1 -ts 1000 500 -co TILED=YES")
gdaltest.runexternal(test_cli_utilities.get_gdalwarp_path() + " data/w_jpeg.tiff tmp/testgdalwarp30_1.tif -t_srs EPSG:3785 -co COMPRESS=LZW -wm 500000 --config GDAL_CACHEMAX 1 -ts 1000 500 -co TILED=YES" + te)

# Second run : with -wo OPTIMIZE_SIZE=TRUE
gdaltest.runexternal(test_cli_utilities.get_gdalwarp_path() + " data/w_jpeg.tiff tmp/testgdalwarp30_2.tif -t_srs EPSG:3785 -co COMPRESS=LZW -wm 500000 -wo OPTIMIZE_SIZE=TRUE --config GDAL_CACHEMAX 1 -ts 1000 500 -co TILED=YES")
gdaltest.runexternal(test_cli_utilities.get_gdalwarp_path() + " data/w_jpeg.tiff tmp/testgdalwarp30_2.tif -t_srs EPSG:3785 -co COMPRESS=LZW -wm 500000 -wo OPTIMIZE_SIZE=TRUE --config GDAL_CACHEMAX 1 -ts 1000 500 -co TILED=YES" + te)

# Third run : with -wo STREAMABLE_OUTPUT=TRUE
gdaltest.runexternal(test_cli_utilities.get_gdalwarp_path() + " data/w_jpeg.tiff tmp/testgdalwarp30_3.tif -t_srs EPSG:3785 -co COMPRESS=LZW -wm 500000 -wo STREAMABLE_OUTPUT=TRUE --config GDAL_CACHEMAX 1 -ts 1000 500 -co TILED=YES")
gdaltest.runexternal(test_cli_utilities.get_gdalwarp_path() + " data/w_jpeg.tiff tmp/testgdalwarp30_3.tif -t_srs EPSG:3785 -co COMPRESS=LZW -wm 500000 -wo STREAMABLE_OUTPUT=TRUE --config GDAL_CACHEMAX 1 -ts 1000 500 -co TILED=YES" + te)

file_size1 = os.stat('tmp/testgdalwarp30_1.tif')[stat.ST_SIZE]
file_size2 = os.stat('tmp/testgdalwarp30_2.tif')[stat.ST_SIZE]
Expand Down
121 changes: 63 additions & 58 deletions gdal/frmts/gtiff/cogdriver.cpp
View file
Expand Up @@ -151,11 +151,70 @@ bool COGGetWarpingCharacteristics(GDALDataset* poSrcDS,

CPLStringList aosTO;
aosTO.SetNameValue( "DST_SRS", osTargetSRS );
void* hTransformArg =
GDALCreateGenImgProjTransformer2( poSrcDS, nullptr, aosTO.List() );
void* hTransformArg = nullptr;

OGRSpatialReference oTargetSRS;
oTargetSRS.SetFromUserInput(osTargetSRS);
const char* pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
const int nEPSGCode = pszAuthCode ? atoi(pszAuthCode) : 0;

// Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
// ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
// EPSG:3857.
double adfSrcGeoTransform[6];
std::unique_ptr<GDALDataset> poTmpDS;
if( nEPSGCode == 3857 &&
poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None &&
adfSrcGeoTransform[2] == 0 &&
adfSrcGeoTransform[4] == 0 &&
adfSrcGeoTransform[5] < 0 )
{
const auto poSrcSRS = poSrcDS->GetSpatialRef();
if( poSrcSRS && poSrcSRS->IsGeographic() )
{
double maxLat = adfSrcGeoTransform[3];
double minLat = adfSrcGeoTransform[3] +
poSrcDS->GetRasterYSize() *
adfSrcGeoTransform[5];
// Corresponds to the latitude of below MAX_GM
constexpr double MAX_LAT = 85.0511287798066;
bool bModified = false;
if( maxLat > MAX_LAT )
{
maxLat = MAX_LAT;
bModified = true;
}
if( minLat < -MAX_LAT )
{
minLat = -MAX_LAT;
bModified = true;
}
if( bModified )
{
CPLStringList aosOptions;
aosOptions.AddString("-of");
aosOptions.AddString("VRT");
aosOptions.AddString("-projwin");
aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0]));
aosOptions.AddString(CPLSPrintf("%.18g", maxLat));
aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0] + poSrcDS->GetRasterXSize() * adfSrcGeoTransform[1]));
aosOptions.AddString(CPLSPrintf("%.18g", minLat));
auto psOptions = GDALTranslateOptionsNew(aosOptions.List(), nullptr);
poTmpDS.reset(GDALDataset::FromHandle(
GDALTranslate("", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
GDALTranslateOptionsFree(psOptions);
if( poTmpDS )
{
hTransformArg = GDALCreateGenImgProjTransformer2(
GDALDataset::FromHandle(poTmpDS.get()), nullptr, aosTO.List() );
}
}
}
}
if( hTransformArg == nullptr )
{
return false;
hTransformArg =
GDALCreateGenImgProjTransformer2( poSrcDS, nullptr, aosTO.List() );
}

GDALTransformerInfo* psInfo = static_cast<GDALTransformerInfo*>(hTransformArg);
Expand All @@ -174,61 +233,7 @@ bool COGGetWarpingCharacteristics(GDALDataset* poSrcDS,

GDALDestroyGenImgProjTransformer( hTransformArg );
hTransformArg = nullptr;


// Hack to compensate for GDALSuggestedWarpOutput2() failure when
// reprojection latitude = +/- 90 to EPSG:3857.
double adfSrcGeoTransform[6];
OGRSpatialReference oTargetSRS;
oTargetSRS.SetFromUserInput(osTargetSRS);
const char* pszAuthCode = oTargetSRS.GetAuthorityCode(nullptr);
const int nEPSGCode = pszAuthCode ? atoi(pszAuthCode) : 0;
if( nEPSGCode == 3857 && poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None )
{
const char* pszSrcWKT = poSrcDS->GetProjectionRef();
if( pszSrcWKT != nullptr && pszSrcWKT[0] != '\0' )
{
OGRSpatialReference oSrcSRS;
if( oSrcSRS.SetFromUserInput( pszSrcWKT ) == OGRERR_NONE &&
oSrcSRS.IsGeographic() )
{
const double minLat =
std::min(adfSrcGeoTransform[3],
adfSrcGeoTransform[3] +
poSrcDS->GetRasterYSize() *
adfSrcGeoTransform[5]);
const double maxLat =
std::max(adfSrcGeoTransform[3],
adfSrcGeoTransform[3] +
poSrcDS->GetRasterYSize() *
adfSrcGeoTransform[5]);
double maxNorthing = adfGeoTransform[3];
double minNorthing =
adfGeoTransform[3] + adfGeoTransform[5] * nYSize;
bool bChanged = false;
const double SPHERICAL_RADIUS = 6378137.0;
const double MAX_GM =
SPHERICAL_RADIUS * M_PI; // 20037508.342789244
if( maxLat > 89.9999999 )
{
bChanged = true;
maxNorthing = MAX_GM;
}
if( minLat <= -89.9999999 )
{
bChanged = true;
minNorthing = -MAX_GM;
}
if( bChanged )
{
adfGeoTransform[3] = maxNorthing;
nYSize = int((maxNorthing - minNorthing) / (-adfGeoTransform[5]) + 0.5);
adfExtent[1] = maxNorthing + nYSize * adfGeoTransform[5];
adfExtent[3] = maxNorthing;
}
}
}
}
poTmpDS.reset();

dfMinX = adfExtent[0];
dfMinY = adfExtent[1];
Expand Down
118 changes: 78 additions & 40 deletions gdal/frmts/mbtiles/mbtilesdataset.cpp
View file
Expand Up @@ -37,6 +37,7 @@
#include "cpl_json.h"
#include "cpl_vsil_curl_priv.h"
#include "gpkgmbtilescommon.h"
#include "gdal_utils.h"
#include "gdalwarper.h"
#include "mvtutils.h"

Expand Down Expand Up @@ -3073,8 +3074,67 @@ GDALDataset* MBTilesDataset::CreateCopy( const char *pszFilename,
}

char** papszTO = CSLSetNameValue( nullptr, "DST_SRS", SRS_EPSG_3857 );
void* hTransformArg =

void* hTransformArg = nullptr;

// Hack to compensate for GDALSuggestedWarpOutput2() failure (or not
// ideal suggestion with PROJ 8) when reprojecting latitude = +/- 90 to
// EPSG:3857.
double adfSrcGeoTransform[6] = {0,0,0,0,0,0};
std::unique_ptr<GDALDataset> poTmpDS;
bool bModifiedMaxLat = false;
bool bModifiedMinLat = false;
const auto poSrcSRS = poSrcDS->GetSpatialRef();
if( poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None &&
adfSrcGeoTransform[2] == 0 &&
adfSrcGeoTransform[4] == 0 &&
adfSrcGeoTransform[5] < 0 )
{
if( poSrcSRS && poSrcSRS->IsGeographic() )
{
double maxLat = adfSrcGeoTransform[3];
double minLat = adfSrcGeoTransform[3] +
poSrcDS->GetRasterYSize() *
adfSrcGeoTransform[5];
// Corresponds to the latitude of MAX_GM
constexpr double MAX_LAT = 85.0511287798066;
if( maxLat > MAX_LAT )
{
maxLat = MAX_LAT;
bModifiedMaxLat = true;
}
if( minLat < -MAX_LAT )
{
minLat = -MAX_LAT;
bModifiedMinLat = true;
}
if( bModifiedMaxLat || bModifiedMinLat )
{
CPLStringList aosOptions;
aosOptions.AddString("-of");
aosOptions.AddString("VRT");
aosOptions.AddString("-projwin");
aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0]));
aosOptions.AddString(CPLSPrintf("%.18g", maxLat));
aosOptions.AddString(CPLSPrintf("%.18g", adfSrcGeoTransform[0] + poSrcDS->GetRasterXSize() * adfSrcGeoTransform[1]));
aosOptions.AddString(CPLSPrintf("%.18g", minLat));
auto psOptions = GDALTranslateOptionsNew(aosOptions.List(), nullptr);
poTmpDS.reset(GDALDataset::FromHandle(
GDALTranslate("", GDALDataset::ToHandle(poSrcDS), psOptions, nullptr)));
GDALTranslateOptionsFree(psOptions);
if( poTmpDS )
{
hTransformArg = GDALCreateGenImgProjTransformer2(
GDALDataset::FromHandle(poTmpDS.get()), nullptr, papszTO );
}
}
}
}
if( hTransformArg == nullptr )
{
hTransformArg =
GDALCreateGenImgProjTransformer2( poSrcDS, nullptr, papszTO );
}
if( hTransformArg == nullptr )
{
CSLDestroy(papszTO);
Expand All @@ -3099,49 +3159,27 @@ GDALDataset* MBTilesDataset::CreateCopy( const char *pszFilename,

GDALDestroyGenImgProjTransformer( hTransformArg );
hTransformArg = nullptr;
poTmpDS.reset();

// Hack to compensate for GDALSuggestedWarpOutput2() failure when
// reprojection latitude = +/- 90 to EPSG:3857
double adfSrcGeoTransform[6];
if( poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None )
if( bModifiedMaxLat || bModifiedMinLat )
{
const char* pszSrcWKT = poSrcDS->GetProjectionRef();
if( pszSrcWKT != nullptr && pszSrcWKT[0] != '\0' )
if( bModifiedMaxLat )
{
const double maxNorthing = MAX_GM;
adfGeoTransform[3] = maxNorthing;
adfExtent[3] = maxNorthing;
}
if( bModifiedMinLat )
{
const double minNorthing = -MAX_GM;
adfExtent[1] = minNorthing;
}

if( poSrcSRS && poSrcSRS->IsGeographic() )
{
OGRSpatialReference oSRS;
if( oSRS.SetFromUserInput( pszSrcWKT ) == OGRERR_NONE &&
oSRS.IsGeographic() )
if( adfSrcGeoTransform[0] + poSrcDS->GetRasterXSize() * adfSrcGeoTransform[1] == 180 )
{
const double minLat =
std::min(
adfSrcGeoTransform[3],
adfSrcGeoTransform[3] +
poSrcDS->GetRasterYSize() * adfSrcGeoTransform[5]);
const double maxLat =
std::max(
adfSrcGeoTransform[3],
adfSrcGeoTransform[3] +
poSrcDS->GetRasterYSize() * adfSrcGeoTransform[5]);
double maxNorthing = adfGeoTransform[3];
double minNorthing = adfGeoTransform[3] + adfGeoTransform[5] * nYSize;
bool bChanged = false;
if( maxLat > 89.9999999 )
{
bChanged = true;
maxNorthing = MAX_GM;
}
if( minLat <= -89.9999999 )
{
bChanged = true;
minNorthing = -MAX_GM;
}
if( bChanged )
{
adfGeoTransform[3] = maxNorthing;
nYSize = int((maxNorthing - minNorthing) / (-adfGeoTransform[5]) + 0.5);
adfExtent[1] = maxNorthing + nYSize * adfGeoTransform[5];
adfExtent[3] = maxNorthing;
}
adfExtent[2] = MAX_GM;
}
}
}
Expand Down

0 comments on commit dc38aa6

Please sign in to comment.