diff --git a/autotest/gdrivers/mbtiles.py b/autotest/gdrivers/mbtiles.py index 93adecf9d9f9..4d8521f2d1e0 100755 --- a/autotest/gdrivers/mbtiles.py +++ b/autotest/gdrivers/mbtiles.py @@ -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() diff --git a/autotest/pyscripts/test_gdal2tiles.py b/autotest/pyscripts/test_gdal2tiles.py index 785eca0d7add..2caaadbc6948 100755 --- a/autotest/pyscripts/test_gdal2tiles.py +++ b/autotest/pyscripts/test_gdal2tiles.py @@ -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']: @@ -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') @@ -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']: diff --git a/autotest/utilities/test_gdalwarp.py b/autotest/utilities/test_gdalwarp.py index 3fd920f88489..e0bd1390edfb 100755 --- a/autotest/utilities/test_gdalwarp.py +++ b/autotest/utilities/test_gdalwarp.py @@ -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() @@ -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] diff --git a/gdal/frmts/gtiff/cogdriver.cpp b/gdal/frmts/gtiff/cogdriver.cpp index 2a6796571f99..795772c9bec2 100644 --- a/gdal/frmts/gtiff/cogdriver.cpp +++ b/gdal/frmts/gtiff/cogdriver.cpp @@ -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 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(hTransformArg); @@ -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]; diff --git a/gdal/frmts/mbtiles/mbtilesdataset.cpp b/gdal/frmts/mbtiles/mbtilesdataset.cpp index 248385f74193..9cf93cc51a6c 100644 --- a/gdal/frmts/mbtiles/mbtilesdataset.cpp +++ b/gdal/frmts/mbtiles/mbtilesdataset.cpp @@ -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" @@ -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 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); @@ -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; } } } diff --git a/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp b/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp index b2471a285570..0eb1bb9982d9 100644 --- a/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp +++ b/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp @@ -31,6 +31,7 @@ #include "ogr_p.h" #include "ogr_swq.h" #include "gdalwarper.h" +#include "gdal_utils.h" #include "ogrgeopackageutility.h" #include "ogrsqliteutility.h" #include "vrt/vrtdataset.h" @@ -4786,8 +4787,70 @@ GDALDataset* GDALGeoPackageDataset::CreateCopy( const char *pszFilename, char* pszWKT = nullptr; oSRS.exportToWkt(&pszWKT); char** papszTO = CSLSetNameValue( nullptr, "DST_SRS", pszWKT ); - 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]; + std::unique_ptr poTmpDS; + bool bEPSG3857Adjust = false; + 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 ) + { + bEPSG3857Adjust = true; + hTransformArg = GDALCreateGenImgProjTransformer2( + GDALDataset::FromHandle(poTmpDS.get()), nullptr, papszTO ); + } + } + } + } + if( hTransformArg == nullptr ) + { + hTransformArg = GDALCreateGenImgProjTransformer2( poSrcDS, nullptr, papszTO ); + } + if( hTransformArg == nullptr ) { CPLFree(pszWKT); @@ -4814,54 +4877,33 @@ GDALDataset* GDALGeoPackageDataset::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( nEPSGCode == 3857 && poSrcDS->GetGeoTransform(adfSrcGeoTransform) == CE_None ) + if( bEPSG3857Adjust ) { - 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 = + constexpr double SPHERICAL_RADIUS = 6378137.0; + constexpr 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; - } - } + double maxNorthing = adfGeoTransform[3]; + double minNorthing = + adfGeoTransform[3] + adfGeoTransform[5] * nYSize; + bool bChanged = false; + if( maxNorthing > MAX_GM ) + { + bChanged = true; + maxNorthing = MAX_GM; + } + if( minNorthing < -MAX_GM ) + { + 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; } } diff --git a/gdal/swig/python/scripts/gdal2tiles.py b/gdal/swig/python/scripts/gdal2tiles.py index 8a2d25151855..bd026c2cd47b 100755 --- a/gdal/swig/python/scripts/gdal2tiles.py +++ b/gdal/swig/python/scripts/gdal2tiles.py @@ -908,6 +908,26 @@ def reproject_dataset(from_dataset, from_srs, to_srs, options=None): raise GDALError("from and to SRS must be defined to reproject the dataset") if (from_srs.ExportToProj4() != to_srs.ExportToProj4()) or (from_dataset.GetGCPCount() != 0): + + if from_srs.IsGeographic() and to_srs.GetAuthorityName(None) == 'EPSG' and to_srs.GetAuthorityCode(None) == '3857': + from_gt = from_dataset.GetGeoTransform(can_return_null=True) + if from_gt and from_gt[2] == 0 and from_gt[4] == 0 and from_gt[5] < 0: + maxlat = from_gt[3] + minlat = from_gt[3] + from_dataset.RasterYSize * from_gt[5] + MAX_LAT = 85.0511287798066 + adjustBounds = False + if maxlat > MAX_LAT: + maxlat = MAX_LAT + adjustBounds = True + if minlat < -MAX_LAT: + minlat = -MAX_LAT + adjustBounds = True + if adjustBounds: + ct = osr.CoordinateTransformation(from_srs, to_srs) + west, south = ct.TransformPoint(from_gt[0], minlat)[:2] + east, north = ct.TransformPoint(from_gt[0] + from_dataset.RasterXSize * from_gt[1], maxlat)[:2] + return gdal.Warp("", from_dataset, format='VRT', outputBounds = [west, south, east, north], dstSRS = 'EPSG:3857') + to_dataset = gdal.AutoCreateWarpedVRT(from_dataset, from_srs.ExportToWkt(), to_srs.ExportToWkt())