38 projector->mSrcCRS = mSrcCRS;
39 projector->mDestCRS = mDestCRS;
40 projector->mTransformContext = mTransformContext;
43 projector->mSrcDatumTransform = mSrcDatumTransform;
44 projector->mDestDatumTransform = mDestDatumTransform;
47 projector->mPrecision = mPrecision;
70 int srcDatumTransform,
71 int destDatumTransform )
76 mSrcDatumTransform = srcDatumTransform;
77 mDestDatumTransform = destDatumTransform;
85 mTransformContext = transformContext;
87 mSrcDatumTransform = -1;
88 mDestDatumTransform = -1;
94 : mApproximate(
false )
95 , mInverseCt( inverseCt )
96 , mDestExtent( extent )
105 , mDestRowsPerMatrixRow( 0.0 )
106 , mDestColsPerMatrixCol( 0.0 )
110 , mSqrTolerance( 0.0 )
128 if ( mExtent.isEmpty() )
130 mExtent = provider->
extent();
135 mDestXRes = mDestExtent.
width() / ( mDestCols );
136 mDestYRes = mDestExtent.height() / ( mDestRows );
142 double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
143 mSqrTolerance = myDestRes * myDestRes;
151 mApproximate =
false;
156 mCPCols = mCPRows = 3;
157 for (
int i = 0; i < mCPRows; i++ )
159 QList<QgsPointXY> myRow;
163 mCPMatrix.insert( i, myRow );
165 QList<bool> myLegalRow;
166 myLegalRow.append(
bool(
false ) );
167 myLegalRow.append(
bool(
false ) );
168 myLegalRow.append(
bool(
false ) );
169 mCPLegalMatrix.insert( i, myLegalRow );
171 for (
int i = 0; i < mCPRows; i++ )
173 calcRow( i, inverseCt );
178 bool myColsOK = checkCols( inverseCt );
181 insertRows( inverseCt );
183 bool myRowsOK = checkRows( inverseCt );
186 insertCols( inverseCt );
188 if ( myColsOK && myRowsOK )
195 if ( mCPRows * mCPCols > 0.25 * mDestRows * mDestCols )
199 mApproximate =
false;
207 QgsDebugMsgLevel( QStringLiteral(
"CPMatrix size: mCPRows = %1 mCPCols = %2" ).arg( mCPRows ).arg( mCPCols ), 4 );
208 mDestRowsPerMatrixRow =
static_cast< double >( mDestRows ) / ( mCPRows - 1 );
209 mDestColsPerMatrixCol =
static_cast< double >( mDestCols ) / ( mCPCols - 1 );
219 calcHelper( 0, pHelperTop );
220 calcHelper( 1, pHelperBottom );
226 mSrcYRes = mSrcExtent.height() / mSrcRows;
227 mSrcXRes = mSrcExtent.width() / mSrcCols;
230 ProjectorData::~ProjectorData()
233 delete[] pHelperBottom;
237 void ProjectorData::calcSrcExtent()
248 mSrcExtent =
QgsRectangle( myPoint.
x(), myPoint.
y(), myPoint.
x(), myPoint.
y() );
249 for (
int i = 0; i < mCPRows; i++ )
251 for (
int j = 0; j < mCPCols ; j++ )
253 myPoint = mCPMatrix[i][j];
254 if ( mCPLegalMatrix[i][j] )
256 mSrcExtent.combineExtentWith( myPoint.
x(), myPoint.
y() );
263 mSrcExtent = mSrcExtent.intersect( mExtent );
274 if ( !mExtent.isEmpty() )
276 if ( mMaxSrcXRes > 0 )
279 double col = std::floor( ( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
280 double x = mExtent.xMinimum() + col * mMaxSrcXRes;
281 mSrcExtent.setXMinimum( x );
283 col = std::ceil( ( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
284 x = mExtent.xMinimum() + col * mMaxSrcXRes;
285 mSrcExtent.setXMaximum( x );
287 if ( mMaxSrcYRes > 0 )
289 double row = std::floor( ( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
290 double y = mExtent.yMaximum() - row * mMaxSrcYRes;
291 mSrcExtent.setYMaximum( y );
293 row = std::ceil( ( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
294 y = mExtent.yMaximum() - row * mMaxSrcYRes;
295 mSrcExtent.setYMinimum( y );
301 QString ProjectorData::cpToString()
304 for (
int i = 0; i < mCPRows; i++ )
308 for (
int j = 0; j < mCPCols; j++ )
311 myString += QLatin1String(
" " );
313 if ( mCPLegalMatrix[i][j] )
319 myString += QLatin1String(
"(-,-)" );
326 void ProjectorData::calcSrcRowsCols()
334 double myMinSize = std::numeric_limits<double>::max();
339 double myDestColsPerMatrixCell =
static_cast< double >( mDestCols ) / mCPCols;
340 double myDestRowsPerMatrixCell =
static_cast< double >( mDestRows ) / mCPRows;
341 QgsDebugMsgLevel( QStringLiteral(
"myDestColsPerMatrixCell = %1 myDestRowsPerMatrixCell = %2" ).arg( myDestColsPerMatrixCell ).arg( myDestRowsPerMatrixCell ), 4 );
342 for (
int i = 0; i < mCPRows - 1; i++ )
344 for (
int j = 0; j < mCPCols - 1; j++ )
349 if ( mCPLegalMatrix[i][j] && mCPLegalMatrix[i][j + 1] && mCPLegalMatrix[i + 1][j] )
351 double mySize = std::sqrt( myPointA.
sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
352 if ( mySize < myMinSize )
355 mySize = std::sqrt( myPointA.
sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
356 if ( mySize < myMinSize )
367 int srcXSize, srcYSize;
370 double srcXRes = srcExtent.
width() / srcXSize;
371 double srcYRes = srcExtent.
height() / srcYSize;
372 myMinSize = std::min( srcXRes, srcYRes );
376 QgsDebugMsg( QStringLiteral(
"Cannot get src extent/size" ) );
385 QgsDebugMsgLevel( QStringLiteral(
"mMaxSrcXRes = %1 mMaxSrcYRes = %2" ).arg( mMaxSrcXRes ).arg( mMaxSrcYRes ), 4 );
387 double myMinXSize = mMaxSrcXRes > myMinSize ? mMaxSrcXRes : myMinSize;
388 double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
389 QgsDebugMsgLevel( QStringLiteral(
"myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ), 4 );
390 QgsDebugMsgLevel( QStringLiteral(
"mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ), 4 );
393 mSrcRows =
static_cast< int >( std::round( mSrcExtent.height() / myMinYSize ) );
394 mSrcCols =
static_cast< int >( std::round( mSrcExtent.width() / myMinXSize ) );
396 QgsDebugMsgLevel( QStringLiteral(
"mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ), 4 );
400 inline void ProjectorData::destPointOnCPMatrix(
int row,
int col,
double *theX,
double *theY )
402 *theX = mDestExtent.xMinimum() + col * mDestExtent.width() / ( mCPCols - 1 );
403 *theY = mDestExtent.yMaximum() - row * mDestExtent.height() / ( mCPRows - 1 );
406 inline int ProjectorData::matrixRow(
int destRow )
408 return static_cast< int >( std::floor( ( destRow + 0.5 ) / mDestRowsPerMatrixRow ) );
410 inline int ProjectorData::matrixCol(
int destCol )
412 return static_cast< int >( std::floor( ( destCol + 0.5 ) / mDestColsPerMatrixCol ) );
415 void ProjectorData::calcHelper(
int matrixRow,
QgsPointXY *points )
418 for (
int myDestCol = 0; myDestCol < mDestCols; myDestCol++ )
420 double myDestX = mDestExtent.xMinimum() + ( myDestCol + 0.5 ) * mDestXRes;
422 int myMatrixCol = matrixCol( myDestCol );
424 double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
426 destPointOnCPMatrix( matrixRow, myMatrixCol, &myDestXMin, &myDestYMin );
427 destPointOnCPMatrix( matrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
429 double xfrac = ( myDestX - myDestXMin ) / ( myDestXMax - myDestXMin );
431 QgsPointXY &mySrcPoint0 = mCPMatrix[matrixRow][myMatrixCol];
432 QgsPointXY &mySrcPoint1 = mCPMatrix[matrixRow][myMatrixCol + 1];
433 double s = mySrcPoint0.
x() + ( mySrcPoint1.
x() - mySrcPoint0.
x() ) * xfrac;
434 double t = mySrcPoint0.
y() + ( mySrcPoint1.
y() - mySrcPoint0.
y() ) * xfrac;
436 points[myDestCol].
setX( s );
437 points[myDestCol].
setY( t );
441 void ProjectorData::nextHelper()
446 pHelperTop = pHelperBottom;
448 calcHelper( mHelperTopRow + 2, pHelperBottom );
452 bool ProjectorData::srcRowCol(
int destRow,
int destCol,
int *srcRow,
int *srcCol )
456 return approximateSrcRowCol( destRow, destCol, srcRow, srcCol );
460 return preciseSrcRowCol( destRow, destCol, srcRow, srcCol );
464 bool ProjectorData::preciseSrcRowCol(
int destRow,
int destCol,
int *srcRow,
int *srcCol )
466 #if 0 // too slow, even if we only run it on debug builds! 468 QgsDebugMsgLevel( QStringLiteral(
"theDestRow = %1 mDestExtent.yMaximum() = %2 mDestYRes = %3" ).arg( destRow ).arg( mDestExtent.yMaximum() ).arg( mDestYRes ), 5 );
472 double x = mDestExtent.xMinimum() + ( destCol + 0.5 ) * mDestXRes;
473 double y = mDestExtent.yMaximum() - ( destRow + 0.5 ) * mDestYRes;
480 if ( mInverseCt.isValid() )
484 mInverseCt.transformInPlace( x, y, z );
496 if ( !mExtent.contains(
QgsPointXY( x, y ) ) )
501 *srcRow =
static_cast< int >( std::floor( ( mSrcExtent.yMaximum() - y ) / mSrcYRes ) );
502 *srcCol =
static_cast< int >( std::floor( ( x - mSrcExtent.xMinimum() ) / mSrcXRes ) );
504 QgsDebugMsgLevel( QStringLiteral(
"mSrcExtent.yMinimum() = %1 mSrcExtent.yMaximum() = %2 mSrcYRes = %3" ).arg( mSrcExtent.yMinimum() ).arg( mSrcExtent.yMaximum() ).arg( mSrcYRes ), 5 );
505 QgsDebugMsgLevel( QStringLiteral(
"theSrcRow = %1 srcCol = %2" ).arg( *srcRow ).arg( *srcCol ), 5 );
512 if ( *srcRow >= mSrcRows )
return false;
513 if ( *srcRow < 0 )
return false;
514 if ( *srcCol >= mSrcCols )
return false;
515 if ( *srcCol < 0 )
return false;
520 bool ProjectorData::approximateSrcRowCol(
int destRow,
int destCol,
int *srcRow,
int *srcCol )
522 int myMatrixRow = matrixRow( destRow );
523 int myMatrixCol = matrixCol( destCol );
525 if ( myMatrixRow > mHelperTopRow )
531 double myDestY = mDestExtent.yMaximum() - ( destRow + 0.5 ) * mDestYRes;
535 double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
537 destPointOnCPMatrix( myMatrixRow + 1, myMatrixCol, &myDestXMin, &myDestYMin );
538 destPointOnCPMatrix( myMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
540 double yfrac = ( myDestY - myDestYMin ) / ( myDestYMax - myDestYMin );
549 double tx = myTop.
x();
550 double ty = myTop.
y();
551 double bx = myBot.
x();
552 double by = myBot.
y();
553 double mySrcX = bx + ( tx - bx ) * yfrac;
554 double mySrcY = by + ( ty - by ) * yfrac;
556 if ( !mExtent.contains(
QgsPointXY( mySrcX, mySrcY ) ) )
563 *srcRow =
static_cast< int >( std::floor( ( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes ) );
564 *srcCol =
static_cast< int >( std::floor( ( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes ) );
569 if ( *srcRow >= mSrcRows )
return false;
570 if ( *srcRow < 0 )
return false;
571 if ( *srcCol >= mSrcCols )
return false;
572 if ( *srcCol < 0 )
return false;
579 for (
int r = 0; r < mCPRows - 1; r++ )
581 QList<QgsPointXY> myRow;
582 QList<bool> myLegalRow;
583 myRow.reserve( mCPCols );
584 myLegalRow.reserve( mCPCols );
585 for (
int c = 0;
c < mCPCols; ++
c )
588 myLegalRow.append(
false );
590 QgsDebugMsgLevel( QStringLiteral(
"insert new row at %1" ).arg( 1 + r * 2 ), 3 );
591 mCPMatrix.insert( 1 + r * 2, myRow );
592 mCPLegalMatrix.insert( 1 + r * 2, myLegalRow );
594 mCPRows += mCPRows - 1;
595 for (
int r = 1; r < mCPRows - 1; r += 2 )
603 for (
int r = 0; r < mCPRows; r++ )
605 for (
int c = 0;
c < mCPCols - 1;
c++ )
608 mCPLegalMatrix[r].insert( 1 +
c * 2,
false );
611 mCPCols += mCPCols - 1;
612 for (
int c = 1;
c < mCPCols - 1;
c += 2 )
621 double myDestX, myDestY;
622 destPointOnCPMatrix( row, col, &myDestX, &myDestY );
628 mCPMatrix[row][col] = ct.
transform( myDestPoint );
629 mCPLegalMatrix[row][col] =
true;
633 mCPLegalMatrix[row][col] =
false;
640 mCPLegalMatrix[row][col] =
false;
647 for (
int i = 0; i < mCPCols; i++ )
649 calcCP( row, i, ct );
658 for (
int i = 0; i < mCPRows; i++ )
660 calcCP( i, col, ct );
673 for (
int c = 0;
c < mCPCols;
c++ )
675 for (
int r = 1; r < mCPRows - 1; r += 2 )
677 double myDestX, myDestY;
678 destPointOnCPMatrix( r,
c, &myDestX, &myDestY );
685 QgsPointXY mySrcApprox( ( mySrcPoint1.
x() + mySrcPoint3.
x() ) / 2, ( mySrcPoint1.
y() + mySrcPoint3.
y() ) / 2 );
686 if ( !mCPLegalMatrix[r - 1][
c] || !mCPLegalMatrix[r][
c] || !mCPLegalMatrix[r + 1][
c] )
694 double mySqrDist = myDestApprox.
sqrDist( myDestPoint );
695 if ( mySqrDist > mSqrTolerance )
718 for (
int r = 0; r < mCPRows; r++ )
720 for (
int c = 1;
c < mCPCols - 1;
c += 2 )
722 double myDestX, myDestY;
723 destPointOnCPMatrix( r,
c, &myDestX, &myDestY );
730 QgsPointXY mySrcApprox( ( mySrcPoint1.
x() + mySrcPoint3.
x() ) / 2, ( mySrcPoint1.
y() + mySrcPoint3.
y() ) / 2 );
731 if ( !mCPLegalMatrix[r][
c - 1] || !mCPLegalMatrix[r][
c] || !mCPLegalMatrix[r][
c + 1] )
739 double mySqrDist = myDestApprox.
sqrDist( myDestPoint );
740 if ( mySqrDist > mSqrTolerance )
764 return tr(
"Approximate" );
766 return tr(
"Exact" );
768 return QStringLiteral(
"Unknown" );
774 QgsDebugMsgLevel( QStringLiteral(
"width = %1 height = %2" ).arg( width ).arg( height ), 4 );
784 if ( ! mSrcCRS.
isValid() || ! mDestCRS.
isValid() || mSrcCRS == mDestCRS )
787 return mInput->
block( bandNo, extent, width, height, feedback );
795 ProjectorData pd( extent, width, height,
mInput, inverseCt, mPrecision, feedback );
800 QgsDebugMsgLevel( QStringLiteral(
"srcExtent:\n%1" ).arg( pd.srcExtent().toString() ), 4 );
801 QgsDebugMsgLevel( QStringLiteral(
"srcCols = %1 srcRows = %2" ).arg( pd.srcCols() ).arg( pd.srcRows() ), 4 );
804 if ( pd.srcRows() <= 0 || pd.srcCols() <= 0 )
810 std::unique_ptr< QgsRasterBlock > inputBlock(
mInput->
block( bandNo, pd.srcExtent(), pd.srcCols(), pd.srcRows(), feedback ) );
811 if ( !inputBlock || inputBlock->isEmpty() )
813 QgsDebugMsg( QStringLiteral(
"No raster data!" ) );
819 std::unique_ptr< QgsRasterBlock > outputBlock(
new QgsRasterBlock( inputBlock->dataType(), width, height ) );
820 if ( inputBlock->hasNoDataValue() )
822 outputBlock->setNoDataValue( inputBlock->noDataValue() );
824 if ( !outputBlock->isValid() )
826 QgsDebugMsg( QStringLiteral(
"Cannot create block" ) );
827 return outputBlock.release();
831 outputBlock->setIsNoData();
845 outputBlock->setIsNoData();
848 for (
int i = 0; i < height; ++i )
852 for (
int j = 0; j < width; ++j )
854 bool inside = pd.srcRowCol( i, j, &srcRow, &srcCol );
855 if ( !inside )
continue;
857 qgssize srcIndex =
static_cast< qgssize >( srcRow * pd.srcCols() + srcCol );
860 if ( doNoData && inputBlock->isNoData( srcRow, srcCol ) )
862 outputBlock->setIsNoData( i, j );
867 char *srcBits = inputBlock->bits( srcIndex );
868 char *destBits = outputBlock->bits( destIndex );
879 memcpy( destBits, srcBits, pixelSize );
880 outputBlock->setIsData( i, j );
884 return outputBlock.release();
888 QgsRectangle &destExtent,
int &destXSize,
int &destYSize )
890 if ( srcExtent.
isEmpty() || srcXSize <= 0 || srcYSize <= 0 )
900 return extentSize( ct, srcExtent, srcXSize, srcYSize, destExtent, destXSize, destYSize );
904 const QgsRectangle &srcExtent,
int srcXSize,
int srcYSize,
905 QgsRectangle &destExtent,
int &destXSize,
int &destYSize )
907 if ( srcExtent.
isEmpty() || srcXSize <= 0 || srcYSize <= 0 )
916 double srcXStep = srcExtent.
width() / 3;
917 double srcYStep = srcExtent.
height() / 3;
918 double srcXRes = srcExtent.
width() / srcXSize;
919 double srcYRes = srcExtent.
height() / srcYSize;
920 double destXRes = std::numeric_limits<double>::max();
921 double destYRes = std::numeric_limits<double>::max();
923 for (
int i = 0; i < 3; i++ )
925 double x = srcExtent.
xMinimum() + i * srcXStep;
926 for (
int j = 0; j < 3; j++ )
928 double y = srcExtent.
yMinimum() + j * srcYStep;
929 QgsRectangle srcRectangle( x - srcXRes / 2, y - srcYRes / 2, x + srcXRes / 2, y + srcYRes / 2 );
933 if ( destRectangle.
width() > 0 )
935 destXRes = std::min( destXRes, destRectangle.
width() );
937 if ( destRectangle.
height() > 0 )
939 destYRes = std::min( destYRes, destRectangle.
height() );
948 destXSize = std::max( 1, static_cast< int >( destExtent.
width() / destYRes ) );
949 destYSize = std::max( 1, static_cast< int >( destExtent.
height() / destYRes ) );
virtual int bandCount() const =0
Gets number of bands.
A rectangle specified with double values.
Precision precision() const
Approximate (default), fast but possibly inaccurate.
virtual QgsRectangle extent() const
Gets the extent of the interface.
A class to represent a 2D point.
QgsRasterProjector * clone() const override
Clone itself, create deep copy.
static bool typeIsNumeric(Qgis::DataType type)
Returns true if data type is numeric.
virtual QgsRasterInterface * input() const
Current input.
#define Q_NOWARN_DEPRECATED_PUSH
DataType
Raster data types.
virtual int ySize() const
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
virtual Qgis::DataType dataType(int bandNo) const =0
Returns data type for the band specified by number.
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
Qgis::DataType dataType(int bandNo) const override
Returns data type for the band specified by number.
#define QgsDebugMsgLevel(str, level)
bool isEmpty() const
Returns true if the rectangle is empty.
virtual QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr)=0
Read block of data using given extent and size.
virtual int capabilities() const
Returns a bitmask containing the supported capabilities.
QgsRectangle extent() const override=0
Returns the extent of the layer.
double width() const
Returns the width of the rectangle.
void setY(double y)
Sets the y value of the point.
QString toString(int precision=16) const
Returns a string representation of form xmin,ymin : xmax,ymax Coordinates will be truncated to the sp...
Unknown or unspecified type.
static int typeSize(int dataType)
Contains information about the context in which a coordinate transform is executed.
Base class for processing filters like renderers, reprojector, resampler etc.
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
void setX(double x)
Sets the x value of the point.
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
static bool extentSize(const QgsCoordinateTransform &ct, const QgsRectangle &srcExtent, int srcXSize, int srcYSize, QgsRectangle &destExtent, int &destXSize, int &destYSize)
Calculate destination extent and size from source extent and size.
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
#define Q_NOWARN_DEPRECATED_POP
Precision
Precision defines if each pixel is reprojected or approximate reprojection based on an approximation ...
bool isCanceled() const
Tells whether the operation has been canceled already.
This class represents a coordinate reference system (CRS).
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Custom exception class for Coordinate Reference System related exceptions.
QgsRasterInterface * mInput
Feedback object tailored for raster block reading.
virtual int xSize() const
Gets raster size.
bool destExtentSize(const QgsRectangle &srcExtent, int srcXSize, int srcYSize, QgsRectangle &destExtent, int &destXSize, int &destYSize)
Calculate destination extent and size from source extent and size.
static QString precisionLabel(Precision precision)
virtual const QgsRasterInterface * sourceInput() const
Gets source / raw input, the first in pipe, usually provider.
double height() const
Returns the height of the rectangle.
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.
int bandCount() const override
Gets number of bands.
Base class for raster data providers.
bool isValid() const
Returns whether this CRS is correctly initialized and usable.