Skip to content

Commit

Permalink
QgsRasterProjection: avoid excessive resolution when reprojecting fro…
Browse files Browse the repository at this point in the history 
…m EPSG:4326 to EPSG:3857, and misalignment issues

Fixes #34518
  • Loading branch information
@rouault @github-actions
rouault authored and github-actions[bot] committed Jan 24, 2022
1 parent 2c71239 commit 9585748
Showing 1 changed file with 65 additions and 9 deletions.
74 changes: 65 additions & 9 deletions src/core/raster/qgsrasterprojector.cpp
View file
Expand Up @@ -159,6 +159,16 @@ ProjectorData::ProjectorData( const QgsRectangle &extent, int width, int height,
// Always try to calculate mCPMatrix, it is used in calcSrcExtent() for both Approximate and Exact
// Initialize the matrix by corners and middle points
mCPCols = mCPRows = 3;

// For WebMercator to geographic, if the bounding box is symmetric in Y,
// and we only use 3 sample points, we would believe that there is perfect
// linear approximation, resulting in wrong reprojection
// (see https://github.com/qgis/QGIS/issues/34518).
// But if we use 5 sample points, we will detect the non-linearity and will
// refine the CPMatrix appropriately.
if ( std::fabs( -mDestExtent.yMinimum() - mDestExtent.yMaximum() ) / height < 0.5 * mDestYRes )
mCPRows = 5;

for ( int i = 0; i < mCPRows; i++ )
{
QList<QgsPointXY> myRow;
Expand Down Expand Up @@ -340,6 +350,8 @@ void ProjectorData::calcSrcRowsCols()

if ( mApproximate )
{
double myMaxSize = 0;

// For now, we take cell sizes projected to source but not to source axes
const double myDestColsPerMatrixCell = static_cast< double >( mDestCols ) / mCPCols;
const double myDestRowsPerMatrixCell = static_cast< double >( mDestRows ) / mCPRows;
Expand All @@ -356,13 +368,22 @@ void ProjectorData::calcSrcRowsCols()
double mySize = std::sqrt( myPointA.sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
if ( mySize < myMinSize )
myMinSize = mySize;
if ( mySize > myMaxSize )
myMaxSize = mySize;

mySize = std::sqrt( myPointA.sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
if ( mySize < myMinSize )
myMinSize = mySize;
if ( mySize > myMaxSize )
myMaxSize = mySize;
}
}
}
// Limit resolution to 1/10th of the maximum resolution to avoid issues
// with for example WebMercator at high northings that result in very small
// latitude differences.
if ( myMinSize < 0.1 * myMaxSize )
myMinSize = 0.1 * myMaxSize;
}
else
{
Expand Down Expand Up @@ -395,8 +416,19 @@ void ProjectorData::calcSrcRowsCols()
QgsDebugMsgLevel( QStringLiteral( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ), 4 );

// we have to round to keep alignment set in calcSrcExtent
mSrcRows = static_cast< int >( std::round( mSrcExtent.height() / myMinYSize ) );
mSrcCols = static_cast< int >( std::round( mSrcExtent.width() / myMinXSize ) );
// Limit to 10x the source dimensions to avoid excessive memory allocation
// and processing time.
double dblSrcRows = mSrcExtent.height() / myMinYSize;
if ( dblSrcRows > mDestRows * 10 )
mSrcRows = mDestRows * 10;
else
mSrcRows = static_cast< int >( std::round( dblSrcRows ) );

double dblSrcCols = mSrcExtent.width() / myMinXSize;
if ( dblSrcCols > mDestCols * 10 )
mSrcCols = mDestCols * 10;
else
mSrcCols = static_cast< int >( std::round( dblSrcCols ) );

QgsDebugMsgLevel( QStringLiteral( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ), 4 );
}
Expand Down Expand Up @@ -859,7 +891,7 @@ QgsRasterBlock *QgsRasterProjector::block( int bandNo, QgsRectangle const &exte
const bool inside = pd.srcRowCol( i, j, &srcRow, &srcCol );
if ( !inside ) continue; // we have everything set to no data

const qgssize srcIndex = static_cast< qgssize >( srcRow * pd.srcCols() + srcCol );
const qgssize srcIndex = static_cast< qgssize >( srcRow ) * pd.srcCols() + srcCol;

// isNoData() may be slow so we check doNoData first
if ( doNoData && inputBlock->isNoData( srcRow, srcCol ) )
Expand All @@ -868,7 +900,7 @@ QgsRasterBlock *QgsRasterProjector::block( int bandNo, QgsRectangle const &exte
continue;
}

const qgssize destIndex = static_cast< qgssize >( i * width + j );
const qgssize destIndex = static_cast< qgssize >( i ) * width + j;
char *srcBits = inputBlock->bits( srcIndex );
char *destBits = outputBlock->bits( destIndex );
if ( !srcBits )
Expand Down Expand Up @@ -918,17 +950,20 @@ bool QgsRasterProjector::extentSize( const QgsCoordinateTransform &ct,

// We reproject pixel rectangle from 9 points matrix of source extent, of course, it gives
// bigger xRes,yRes than reprojected edges (envelope)
const double srcXStep = srcExtent.width() / 3;
const double srcYStep = srcExtent.height() / 3;
constexpr int steps = 3;
const double srcXStep = srcExtent.width() / steps;
const double srcYStep = srcExtent.height() / steps;
const double srcXRes = srcExtent.width() / srcXSize;
const double srcYRes = srcExtent.height() / srcYSize;
double destXRes = std::numeric_limits<double>::max();
double destYRes = std::numeric_limits<double>::max();
double maxXRes = 0;
double maxYRes = 0;

for ( int i = 0; i < 3; i++ )
for ( int i = 0; i < steps; i++ )
{
const double x = srcExtent.xMinimum() + i * srcXStep;
for ( int j = 0; j < 3; j++ )
for ( int j = 0; j < steps; j++ )
{
const double y = srcExtent.yMinimum() + j * srcYStep;
const QgsRectangle srcRectangle( x - srcXRes / 2, y - srcYRes / 2, x + srcXRes / 2, y + srcYRes / 2 );
Expand All @@ -938,10 +973,14 @@ bool QgsRasterProjector::extentSize( const QgsCoordinateTransform &ct,
if ( destRectangle.width() > 0 )
{
destXRes = std::min( destXRes, destRectangle.width() );
if ( destRectangle.width() > maxXRes )
maxXRes = destRectangle.width();
}
if ( destRectangle.height() > 0 )
{
destYRes = std::min( destYRes, destRectangle.height() );
if ( destRectangle.height() > maxYRes )
maxYRes = destRectangle.height();
}
}
catch ( QgsCsException & )
Expand All @@ -950,7 +989,24 @@ bool QgsRasterProjector::extentSize( const QgsCoordinateTransform &ct,
}
}
}
destXSize = std::max( 1, static_cast< int >( destExtent.width() / destYRes ) );

// Limit resolution to 1/10th of the maximum resolution to avoid issues
// with for example WebMercator at high northings that result in very small
// latitude differences.
if ( destXRes < 0.1 * maxXRes )
{
destXRes = 0.1 * maxXRes;
}
if ( destYRes < 0.1 * maxYRes )
{
destYRes = 0.1 * maxYRes;
}
if ( destXRes == 0 || destExtent.width() / destXRes > std::numeric_limits<int>::max() )
return false;
if ( destYRes == 0 || destExtent.height() / destYRes > std::numeric_limits<int>::max() )
return false;

destXSize = std::max( 1, static_cast< int >( destExtent.width() / destXRes ) );
destYSize = std::max( 1, static_cast< int >( destExtent.height() / destYRes ) );

return true;
Expand Down

0 comments on commit 9585748

Please sign in to comment.