QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
qgsrasterprojector.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsrasterprojector.cpp - Raster projector
3 --------------------------------------
4 Date : Jan 16, 2011
5 Copyright : (C) 2005 by Radim Blazek
6 email : radim dot blazek at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17#include <algorithm>
18
20#include "qgslogger.h"
21#include "qgsrasterprojector.h"
23#include "qgsexception.h"
24
25Q_NOWARN_DEPRECATED_PUSH // because of deprecated members
27 : QgsRasterInterface( nullptr )
28{
29 QgsDebugMsgLevel( QStringLiteral( "Entered" ), 4 );
30}
32
33
35{
36 QgsDebugMsgLevel( QStringLiteral( "Entered" ), 4 );
38 projector->mSrcCRS = mSrcCRS;
39 projector->mDestCRS = mDestCRS;
40 projector->mTransformContext = mTransformContext;
41
43 projector->mSrcDatumTransform = mSrcDatumTransform;
44 projector->mDestDatumTransform = mDestDatumTransform;
46
47 projector->mPrecision = mPrecision;
48 return projector;
49}
50
52{
53 if ( mInput ) return mInput->bandCount();
54
55 return 0;
56}
57
59{
60 if ( mInput ) return mInput->dataType( bandNo );
61
63}
64
65
67
69 const QgsCoordinateReferenceSystem &destCRS,
70 int srcDatumTransform,
71 int destDatumTransform )
72{
73 mSrcCRS = srcCRS;
74 mDestCRS = destCRS;
76 mSrcDatumTransform = srcDatumTransform;
77 mDestDatumTransform = destDatumTransform;
79}
80
82{
83 mSrcCRS = srcCRS;
84 mDestCRS = destCRS;
85 mTransformContext = transformContext;
87 mSrcDatumTransform = -1;
88 mDestDatumTransform = -1;
90}
91
92
93ProjectorData::ProjectorData( const QgsRectangle &extent, int width, int height, QgsRasterInterface *input, const QgsCoordinateTransform &inverseCt, QgsRasterProjector::Precision precision, QgsRasterBlockFeedback *feedback )
94 : mApproximate( false )
95 , mInverseCt( inverseCt )
96 , mDestExtent( extent )
97 , mDestRows( height )
98 , mDestCols( width )
99 , mDestXRes( 0.0 )
100 , mDestYRes( 0.0 )
101 , mSrcRows( 0 )
102 , mSrcCols( 0 )
103 , mSrcXRes( 0.0 )
104 , mSrcYRes( 0.0 )
105 , mDestRowsPerMatrixRow( 0.0 )
106 , mDestColsPerMatrixCol( 0.0 )
107 , mHelperTopRow( 0 )
108 , mCPCols( 0 )
109 , mCPRows( 0 )
110 , mSqrTolerance( 0.0 )
111 , mMaxSrcXRes( 0 )
112 , mMaxSrcYRes( 0 )
113{
114 QgsDebugMsgLevel( QStringLiteral( "Entered" ), 4 );
115
116 // Get max source resolution and extent if possible
117 if ( input )
118 {
119 QgsRasterDataProvider *provider = dynamic_cast<QgsRasterDataProvider *>( input->sourceInput() );
120 if ( provider )
121 {
122 // If provider-side resampling is possible, we will get a much better looking
123 // result by not requesting at the maximum resolution and then doing nearest
124 // resampling here. A real fix would be to do resampling during reprojection
125 // however.
127 ( provider->capabilities() & QgsRasterDataProvider::Size ) )
128 {
129 mMaxSrcXRes = provider->extent().width() / provider->xSize();
130 mMaxSrcYRes = provider->extent().height() / provider->ySize();
131 }
132 // Get source extent
133 if ( mExtent.isEmpty() )
134 {
135 mExtent = provider->extent();
136 }
137 }
138 }
139
140 mDestXRes = mDestExtent.width() / ( mDestCols );
141 mDestYRes = mDestExtent.height() / ( mDestRows );
142
143 // Calculate tolerance
144 // TODO: Think it over better
145 // Note: we are checking on matrix each even point, that means that the real error
146 // in that moment is approximately half size
147 const double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
148 mSqrTolerance = myDestRes * myDestRes;
149
151 {
152 mApproximate = true;
153 }
154 else
155 {
156 mApproximate = false;
157 }
158
159 // Always try to calculate mCPMatrix, it is used in calcSrcExtent() for both Approximate and Exact
160 // Initialize the matrix by corners and middle points
161 mCPCols = mCPRows = 3;
162
163 // For WebMercator to geographic, if the bounding box is symmetric in Y,
164 // and we only use 3 sample points, we would believe that there is perfect
165 // linear approximation, resulting in wrong reprojection
166 // (see https://github.com/qgis/QGIS/issues/34518).
167 // But if we use 5 sample points, we will detect the non-linearity and will
168 // refine the CPMatrix appropriately.
169 if ( std::fabs( -mDestExtent.yMinimum() - mDestExtent.yMaximum() ) / height < 0.5 * mDestYRes )
170 mCPRows = 5;
171
172 for ( int i = 0; i < mCPRows; i++ )
173 {
174 QList<QgsPointXY> myRow;
175 myRow.append( QgsPointXY() );
176 myRow.append( QgsPointXY() );
177 myRow.append( QgsPointXY() );
178 mCPMatrix.insert( i, myRow );
179 // And the legal points
180 QList<bool> myLegalRow;
181 myLegalRow.append( bool( false ) );
182 myLegalRow.append( bool( false ) );
183 myLegalRow.append( bool( false ) );
184 mCPLegalMatrix.insert( i, myLegalRow );
185 }
186 for ( int i = 0; i < mCPRows; i++ )
187 {
188 calcRow( i, inverseCt );
189 }
190
191 while ( true )
192 {
193 const bool myColsOK = checkCols( inverseCt );
194 if ( !myColsOK )
195 {
196 insertRows( inverseCt );
197 }
198 const bool myRowsOK = checkRows( inverseCt );
199 if ( !myRowsOK )
200 {
201 insertCols( inverseCt );
202 }
203 if ( myColsOK && myRowsOK )
204 {
205 QgsDebugMsgLevel( QStringLiteral( "CP matrix within tolerance" ), 4 );
206 break;
207 }
208 // What is the maximum reasonable size of transformatio matrix?
209 // TODO: consider better when to break - ratio
210 if ( mCPRows * mCPCols > 0.25 * mDestRows * mDestCols )
211 //if ( mCPRows * mCPCols > mDestRows * mDestCols )
212 {
213 QgsDebugMsgLevel( QStringLiteral( "Too large CP matrix" ), 4 );
214 mApproximate = false;
215 break;
216 }
217 if ( feedback && feedback->isCanceled() )
218 {
219 return;
220 }
221 }
222 QgsDebugMsgLevel( QStringLiteral( "CPMatrix size: mCPRows = %1 mCPCols = %2" ).arg( mCPRows ).arg( mCPCols ), 4 );
223 mDestRowsPerMatrixRow = static_cast< double >( mDestRows ) / ( mCPRows - 1 );
224 mDestColsPerMatrixCol = static_cast< double >( mDestCols ) / ( mCPCols - 1 );
225
226#if 0
227 QgsDebugMsgLevel( QStringLiteral( "CPMatrix:" ), 5 );
228 QgsDebugMsgLevel( cpToString(), 5 );
229#endif
230
231 // init helper points
232 pHelperTop = new QgsPointXY[mDestCols];
233 pHelperBottom = new QgsPointXY[mDestCols];
234 calcHelper( 0, pHelperTop );
235 calcHelper( 1, pHelperBottom );
236 mHelperTopRow = 0;
237
238 // Calculate source dimensions
239 calcSrcExtent();
240 calcSrcRowsCols();
241 mSrcYRes = mSrcExtent.height() / mSrcRows;
242 mSrcXRes = mSrcExtent.width() / mSrcCols;
243}
244
245ProjectorData::~ProjectorData()
246{
247 delete[] pHelperTop;
248 delete[] pHelperBottom;
249}
250
251
252void ProjectorData::calcSrcExtent()
253{
254 /* Run around the mCPMatrix and find source extent */
255 // Attention, source limits are not necessarily on destination edges, e.g.
256 // for destination EPSG:32661 Polar Stereographic and source EPSG:4326,
257 // the maximum y may be in the middle of destination extent
258 // TODO: How to find extent exactly and quickly?
259 // For now, we run through all matrix
260 // mCPMatrix is used for both Approximate and Exact because QgsCoordinateTransform::transformBoundingBox()
261 // is not precise enough, see #13665
262 QgsPointXY myPoint = mCPMatrix[0][0];
263 mSrcExtent = QgsRectangle( myPoint.x(), myPoint.y(), myPoint.x(), myPoint.y() );
264 for ( int i = 0; i < mCPRows; i++ )
265 {
266 for ( int j = 0; j < mCPCols ; j++ )
267 {
268 myPoint = mCPMatrix[i][j];
269 if ( mCPLegalMatrix[i][j] )
270 {
271 mSrcExtent.combineExtentWith( myPoint.x(), myPoint.y() );
272 }
273 }
274 }
275 // Expand a bit to avoid possible approx coords falling out because of representation error?
276
277 // Combine with maximum source extent
278 mSrcExtent = mSrcExtent.intersect( mExtent );
279
280 // If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
281 // align extent to src resolution to avoid jumping of reprojected pixels
282 // when shifting resampled grid.
283 // Important especially if we are over mMaxSrcXRes, mMaxSrcYRes limits
284 // Note however, that preceding filters (like resampler) may read data
285 // on different resolution.
286
287 QgsDebugMsgLevel( "mSrcExtent = " + mSrcExtent.toString(), 4 );
288 QgsDebugMsgLevel( "mExtent = " + mExtent.toString(), 4 );
289 if ( !mExtent.isEmpty() )
290 {
291 if ( mMaxSrcXRes > 0 )
292 {
293 // with floor/ceil it should work correctly also for mSrcExtent.xMinimum() < mExtent.xMinimum()
294 double col = std::floor( ( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
295 double x = mExtent.xMinimum() + col * mMaxSrcXRes;
296 mSrcExtent.setXMinimum( x );
297
298 col = std::ceil( ( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
299 x = mExtent.xMinimum() + col * mMaxSrcXRes;
300 mSrcExtent.setXMaximum( x );
301 }
302 if ( mMaxSrcYRes > 0 )
303 {
304 double row = std::floor( ( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
305 double y = mExtent.yMaximum() - row * mMaxSrcYRes;
306 mSrcExtent.setYMaximum( y );
307
308 row = std::ceil( ( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
309 y = mExtent.yMaximum() - row * mMaxSrcYRes;
310 mSrcExtent.setYMinimum( y );
311 }
312 }
313 QgsDebugMsgLevel( "mSrcExtent = " + mSrcExtent.toString(), 4 );
314}
315
316QString ProjectorData::cpToString() const
317{
318 QString myString;
319 for ( int i = 0; i < mCPRows; i++ )
320 {
321 if ( i > 0 )
322 myString += '\n';
323 for ( int j = 0; j < mCPCols; j++ )
324 {
325 if ( j > 0 )
326 myString += QLatin1String( " " );
327 const QgsPointXY myPoint = mCPMatrix[i][j];
328 if ( mCPLegalMatrix[i][j] )
329 {
330 myString += myPoint.toString();
331 }
332 else
333 {
334 myString += QLatin1String( "(-,-)" );
335 }
336 }
337 }
338 return myString;
339}
340
341void ProjectorData::calcSrcRowsCols()
342{
343 // Wee need to calculate minimum cell size in the source
344 // TODO: Think it over better, what is the right source resolution?
345 // Taking distances between cell centers projected to source along source
346 // axis would result in very high resolution
347 // TODO: different resolution for rows and cols ?
348
349 double myMinSize = std::numeric_limits<double>::max();
350
351 if ( mApproximate )
352 {
353 double myMaxSize = 0;
354
355 // For now, we take cell sizes projected to source but not to source axes
356 const double myDestColsPerMatrixCell = static_cast< double >( mDestCols ) / mCPCols;
357 const double myDestRowsPerMatrixCell = static_cast< double >( mDestRows ) / mCPRows;
358 QgsDebugMsgLevel( QStringLiteral( "myDestColsPerMatrixCell = %1 myDestRowsPerMatrixCell = %2" ).arg( myDestColsPerMatrixCell ).arg( myDestRowsPerMatrixCell ), 4 );
359 for ( int i = 0; i < mCPRows - 1; i++ )
360 {
361 for ( int j = 0; j < mCPCols - 1; j++ )
362 {
363 const QgsPointXY myPointA = mCPMatrix[i][j];
364 const QgsPointXY myPointB = mCPMatrix[i][j + 1];
365 const QgsPointXY myPointC = mCPMatrix[i + 1][j];
366 if ( mCPLegalMatrix[i][j] && mCPLegalMatrix[i][j + 1] && mCPLegalMatrix[i + 1][j] )
367 {
368 double mySize = std::sqrt( myPointA.sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
369 if ( mySize < myMinSize )
370 myMinSize = mySize;
371 if ( mySize > myMaxSize )
372 myMaxSize = mySize;
373
374 mySize = std::sqrt( myPointA.sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
375 if ( mySize < myMinSize )
376 myMinSize = mySize;
377 if ( mySize > myMaxSize )
378 myMaxSize = mySize;
379 }
380 }
381 }
382 // Limit resolution to 1/10th of the maximum resolution to avoid issues
383 // with for example WebMercator at high northings that result in very small
384 // latitude differences.
385 if ( myMinSize < 0.1 * myMaxSize )
386 myMinSize = 0.1 * myMaxSize;
387 }
388 else
389 {
390 // take highest from corners, points in in the middle of corners and center (3 x 3 )
391 //double
392 QgsRectangle srcExtent;
393 int srcXSize, srcYSize;
394 if ( QgsRasterProjector::extentSize( mInverseCt, mDestExtent, mDestCols, mDestRows, srcExtent, srcXSize, srcYSize ) )
395 {
396 const double srcXRes = srcExtent.width() / srcXSize;
397 const double srcYRes = srcExtent.height() / srcYSize;
398 myMinSize = std::min( srcXRes, srcYRes );
399 }
400 else
401 {
402 QgsDebugMsg( QStringLiteral( "Cannot get src extent/size" ) );
403 }
404 }
405
406 // Make it a bit higher resolution
407 // TODO: find the best coefficient, attention, increasing resolution for WMS
408 // is changing WMS content
409 myMinSize *= 0.75;
410
411 QgsDebugMsgLevel( QStringLiteral( "mMaxSrcXRes = %1 mMaxSrcYRes = %2" ).arg( mMaxSrcXRes ).arg( mMaxSrcYRes ), 4 );
412 // mMaxSrcXRes, mMaxSrcYRes may be 0 - no limit (WMS)
413 const double myMinXSize = mMaxSrcXRes > myMinSize ? mMaxSrcXRes : myMinSize;
414 const double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
415 QgsDebugMsgLevel( QStringLiteral( "myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ), 4 );
416 QgsDebugMsgLevel( QStringLiteral( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ), 4 );
417
418 // we have to round to keep alignment set in calcSrcExtent
419 // Limit to 10x the source dimensions to avoid excessive memory allocation
420 // and processing time.
421 double dblSrcRows = mSrcExtent.height() / myMinYSize;
422 if ( dblSrcRows > mDestRows * 10 )
423 mSrcRows = mDestRows * 10;
424 else
425 mSrcRows = static_cast< int >( std::round( dblSrcRows ) );
426
427 double dblSrcCols = mSrcExtent.width() / myMinXSize;
428 if ( dblSrcCols > mDestCols * 10 )
429 mSrcCols = mDestCols * 10;
430 else
431 mSrcCols = static_cast< int >( std::round( dblSrcCols ) );
432
433 QgsDebugMsgLevel( QStringLiteral( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ), 4 );
434}
435
436
437inline void ProjectorData::destPointOnCPMatrix( int row, int col, double *theX, double *theY )
438{
439 *theX = mDestExtent.xMinimum() + col * mDestExtent.width() / ( mCPCols - 1 );
440 *theY = mDestExtent.yMaximum() - row * mDestExtent.height() / ( mCPRows - 1 );
441}
442
443inline int ProjectorData::matrixRow( int destRow )
444{
445 return static_cast< int >( std::floor( ( destRow + 0.5 ) / mDestRowsPerMatrixRow ) );
446}
447inline int ProjectorData::matrixCol( int destCol )
448{
449 return static_cast< int >( std::floor( ( destCol + 0.5 ) / mDestColsPerMatrixCol ) );
450}
451
452void ProjectorData::calcHelper( int matrixRow, QgsPointXY *points )
453{
454 // TODO?: should we also precalc dest cell center coordinates for x and y?
455 for ( int myDestCol = 0; myDestCol < mDestCols; myDestCol++ )
456 {
457 const double myDestX = mDestExtent.xMinimum() + ( myDestCol + 0.5 ) * mDestXRes;
458
459 const int myMatrixCol = matrixCol( myDestCol );
460
461 double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
462
463 destPointOnCPMatrix( matrixRow, myMatrixCol, &myDestXMin, &myDestYMin );
464 destPointOnCPMatrix( matrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
465
466 const double xfrac = ( myDestX - myDestXMin ) / ( myDestXMax - myDestXMin );
467
468 const QgsPointXY &mySrcPoint0 = mCPMatrix[matrixRow][myMatrixCol];
469 const QgsPointXY &mySrcPoint1 = mCPMatrix[matrixRow][myMatrixCol + 1];
470 const double s = mySrcPoint0.x() + ( mySrcPoint1.x() - mySrcPoint0.x() ) * xfrac;
471 const double t = mySrcPoint0.y() + ( mySrcPoint1.y() - mySrcPoint0.y() ) * xfrac;
472
473 points[myDestCol].setX( s );
474 points[myDestCol].setY( t );
475 }
476}
477
478void ProjectorData::nextHelper()
479{
480 // We just switch pHelperTop and pHelperBottom, memory is not lost
481 QgsPointXY *tmp = nullptr;
482 tmp = pHelperTop;
483 pHelperTop = pHelperBottom;
484 pHelperBottom = tmp;
485 calcHelper( mHelperTopRow + 2, pHelperBottom );
486 mHelperTopRow++;
487}
488
489bool ProjectorData::srcRowCol( int destRow, int destCol, int *srcRow, int *srcCol )
490{
491 if ( mApproximate )
492 {
493 return approximateSrcRowCol( destRow, destCol, srcRow, srcCol );
494 }
495 else
496 {
497 return preciseSrcRowCol( destRow, destCol, srcRow, srcCol );
498 }
499}
500
501bool ProjectorData::preciseSrcRowCol( int destRow, int destCol, int *srcRow, int *srcCol )
502{
503#if 0 // too slow, even if we only run it on debug builds!
504 QgsDebugMsgLevel( QStringLiteral( "theDestRow = %1" ).arg( destRow ), 5 );
505 QgsDebugMsgLevel( QStringLiteral( "theDestRow = %1 mDestExtent.yMaximum() = %2 mDestYRes = %3" ).arg( destRow ).arg( mDestExtent.yMaximum() ).arg( mDestYRes ), 5 );
506#endif
507
508 // Get coordinate of center of destination cell
509 double x = mDestExtent.xMinimum() + ( destCol + 0.5 ) * mDestXRes;
510 double y = mDestExtent.yMaximum() - ( destRow + 0.5 ) * mDestYRes;
511 double z = 0;
512
513#if 0
514 QgsDebugMsgLevel( QStringLiteral( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
515#endif
516
517 if ( mInverseCt.isValid() )
518 {
519 try
520 {
521 mInverseCt.transformInPlace( x, y, z );
522 }
523 catch ( QgsCsException & )
524 {
525 return false;
526 }
527 }
528
529#if 0
530 QgsDebugMsgLevel( QStringLiteral( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
531#endif
532
533 if ( !mExtent.contains( x, y ) )
534 {
535 return false;
536 }
537 // Get source row col
538 *srcRow = static_cast< int >( std::floor( ( mSrcExtent.yMaximum() - y ) / mSrcYRes ) );
539 *srcCol = static_cast< int >( std::floor( ( x - mSrcExtent.xMinimum() ) / mSrcXRes ) );
540#if 0
541 QgsDebugMsgLevel( QStringLiteral( "mSrcExtent.yMinimum() = %1 mSrcExtent.yMaximum() = %2 mSrcYRes = %3" ).arg( mSrcExtent.yMinimum() ).arg( mSrcExtent.yMaximum() ).arg( mSrcYRes ), 5 );
542 QgsDebugMsgLevel( QStringLiteral( "theSrcRow = %1 srcCol = %2" ).arg( *srcRow ).arg( *srcCol ), 5 );
543#endif
544
545 // With epsg 32661 (Polar Stereographic) it was happening that *srcCol == mSrcCols
546 // For now silently correct limits to avoid crashes
547 // TODO: review
548 // should not happen
549 if ( *srcRow >= mSrcRows ) return false;
550 if ( *srcRow < 0 ) return false;
551 if ( *srcCol >= mSrcCols ) return false;
552 if ( *srcCol < 0 ) return false;
553
554 return true;
555}
556
557bool ProjectorData::approximateSrcRowCol( int destRow, int destCol, int *srcRow, int *srcCol )
558{
559 const int myMatrixRow = matrixRow( destRow );
560 const int myMatrixCol = matrixCol( destCol );
561
562 if ( myMatrixRow > mHelperTopRow )
563 {
564 // TODO: make it more robust (for random, not sequential reading)
565 nextHelper();
566 }
567
568 const double myDestY = mDestExtent.yMaximum() - ( destRow + 0.5 ) * mDestYRes;
569
570 // See the schema in javax.media.jai.WarpGrid doc (but up side down)
571 // TODO: use some kind of cache of values which can be reused
572 double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
573
574 destPointOnCPMatrix( myMatrixRow + 1, myMatrixCol, &myDestXMin, &myDestYMin );
575 destPointOnCPMatrix( myMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
576
577 const double yfrac = ( myDestY - myDestYMin ) / ( myDestYMax - myDestYMin );
578
579 const QgsPointXY &myTop = pHelperTop[destCol];
580 const QgsPointXY &myBot = pHelperBottom[destCol];
581
582 // Warning: this is very SLOW compared to the following code!:
583 //double mySrcX = myBot.x() + (myTop.x() - myBot.x()) * yfrac;
584 //double mySrcY = myBot.y() + (myTop.y() - myBot.y()) * yfrac;
585
586 const double tx = myTop.x();
587 const double ty = myTop.y();
588 const double bx = myBot.x();
589 const double by = myBot.y();
590 const double mySrcX = bx + ( tx - bx ) * yfrac;
591 const double mySrcY = by + ( ty - by ) * yfrac;
592
593 if ( !mExtent.contains( mySrcX, mySrcY ) )
594 {
595 return false;
596 }
597
598 // TODO: check again cell selection (coor is in the middle)
599
600 *srcRow = static_cast< int >( std::floor( ( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes ) );
601 *srcCol = static_cast< int >( std::floor( ( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes ) );
602
603 // For now silently correct limits to avoid crashes
604 // TODO: review
605 // should not happen
606 if ( *srcRow >= mSrcRows ) return false;
607 if ( *srcRow < 0 ) return false;
608 if ( *srcCol >= mSrcCols ) return false;
609 if ( *srcCol < 0 ) return false;
610
611 return true;
612}
613
614void ProjectorData::insertRows( const QgsCoordinateTransform &ct )
615{
616 for ( int r = 0; r < mCPRows - 1; r++ )
617 {
618 QList<QgsPointXY> myRow;
619 QList<bool> myLegalRow;
620 myRow.reserve( mCPCols );
621 myLegalRow.reserve( mCPCols );
622 for ( int c = 0; c < mCPCols; ++c )
623 {
624 myRow.append( QgsPointXY() );
625 myLegalRow.append( false );
626 }
627 QgsDebugMsgLevel( QStringLiteral( "insert new row at %1" ).arg( 1 + r * 2 ), 3 );
628 mCPMatrix.insert( 1 + r * 2, myRow );
629 mCPLegalMatrix.insert( 1 + r * 2, myLegalRow );
630 }
631 mCPRows += mCPRows - 1;
632 for ( int r = 1; r < mCPRows - 1; r += 2 )
633 {
634 calcRow( r, ct );
635 }
636}
637
638void ProjectorData::insertCols( const QgsCoordinateTransform &ct )
639{
640 for ( int r = 0; r < mCPRows; r++ )
641 {
642 for ( int c = 0; c < mCPCols - 1; c++ )
643 {
644 mCPMatrix[r].insert( 1 + c * 2, QgsPointXY() );
645 mCPLegalMatrix[r].insert( 1 + c * 2, false );
646 }
647 }
648 mCPCols += mCPCols - 1;
649 for ( int c = 1; c < mCPCols - 1; c += 2 )
650 {
651 calcCol( c, ct );
652 }
653
654}
655
656void ProjectorData::calcCP( int row, int col, const QgsCoordinateTransform &ct )
657{
658 double myDestX, myDestY;
659 destPointOnCPMatrix( row, col, &myDestX, &myDestY );
660 const QgsPointXY myDestPoint( myDestX, myDestY );
661 try
662 {
663 if ( ct.isValid() )
664 {
665 mCPMatrix[row][col] = ct.transform( myDestPoint );
666 mCPLegalMatrix[row][col] = true;
667 }
668 else
669 {
670 mCPLegalMatrix[row][col] = false;
671 }
672 }
673 catch ( QgsCsException &e )
674 {
675 Q_UNUSED( e )
676 // Caught an error in transform
677 mCPLegalMatrix[row][col] = false;
678 }
679}
680
681bool ProjectorData::calcRow( int row, const QgsCoordinateTransform &ct )
682{
683 QgsDebugMsgLevel( QStringLiteral( "theRow = %1" ).arg( row ), 3 );
684 for ( int i = 0; i < mCPCols; i++ )
685 {
686 calcCP( row, i, ct );
687 }
688
689 return true;
690}
691
692bool ProjectorData::calcCol( int col, const QgsCoordinateTransform &ct )
693{
694 QgsDebugMsgLevel( QStringLiteral( "theCol = %1" ).arg( col ), 3 );
695 for ( int i = 0; i < mCPRows; i++ )
696 {
697 calcCP( i, col, ct );
698 }
699
700 return true;
701}
702
703bool ProjectorData::checkCols( const QgsCoordinateTransform &ct )
704{
705 if ( !ct.isValid() )
706 {
707 return false;
708 }
709
710 for ( int c = 0; c < mCPCols; c++ )
711 {
712 for ( int r = 1; r < mCPRows - 1; r += 2 )
713 {
714 double myDestX, myDestY;
715 destPointOnCPMatrix( r, c, &myDestX, &myDestY );
716 const QgsPointXY myDestPoint( myDestX, myDestY );
717
718 const QgsPointXY mySrcPoint1 = mCPMatrix[r - 1][c];
719 const QgsPointXY mySrcPoint2 = mCPMatrix[r][c];
720 const QgsPointXY mySrcPoint3 = mCPMatrix[r + 1][c];
721
722 const QgsPointXY mySrcApprox( ( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
723 if ( !mCPLegalMatrix[r - 1][c] || !mCPLegalMatrix[r][c] || !mCPLegalMatrix[r + 1][c] )
724 {
725 // There was an error earlier in transform, just abort
726 return false;
727 }
728 try
729 {
730 const QgsPointXY myDestApprox = ct.transform( mySrcApprox, Qgis::TransformDirection::Reverse );
731 const double mySqrDist = myDestApprox.sqrDist( myDestPoint );
732 if ( mySqrDist > mSqrTolerance )
733 {
734 return false;
735 }
736 }
737 catch ( QgsCsException &e )
738 {
739 Q_UNUSED( e )
740 // Caught an error in transform
741 return false;
742 }
743 }
744 }
745 return true;
746}
747
748bool ProjectorData::checkRows( const QgsCoordinateTransform &ct )
749{
750 if ( !ct.isValid() )
751 {
752 return false;
753 }
754
755 for ( int r = 0; r < mCPRows; r++ )
756 {
757 for ( int c = 1; c < mCPCols - 1; c += 2 )
758 {
759 double myDestX, myDestY;
760 destPointOnCPMatrix( r, c, &myDestX, &myDestY );
761
762 const QgsPointXY myDestPoint( myDestX, myDestY );
763 const QgsPointXY mySrcPoint1 = mCPMatrix[r][c - 1];
764 const QgsPointXY mySrcPoint2 = mCPMatrix[r][c];
765 const QgsPointXY mySrcPoint3 = mCPMatrix[r][c + 1];
766
767 const QgsPointXY mySrcApprox( ( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
768 if ( !mCPLegalMatrix[r][c - 1] || !mCPLegalMatrix[r][c] || !mCPLegalMatrix[r][c + 1] )
769 {
770 // There was an error earlier in transform, just abort
771 return false;
772 }
773 try
774 {
775 const QgsPointXY myDestApprox = ct.transform( mySrcApprox, Qgis::TransformDirection::Reverse );
776 const double mySqrDist = myDestApprox.sqrDist( myDestPoint );
777 if ( mySqrDist > mSqrTolerance )
778 {
779 return false;
780 }
781 }
782 catch ( QgsCsException &e )
783 {
784 Q_UNUSED( e )
785 // Caught an error in transform
786 return false;
787 }
788 }
789 }
790 return true;
791}
792
794
795
797{
798 switch ( precision )
799 {
800 case Approximate:
801 return tr( "Approximate" );
802 case Exact:
803 return tr( "Exact" );
804 }
805 return QStringLiteral( "Unknown" );
806}
807
808QgsRasterBlock *QgsRasterProjector::block( int bandNo, QgsRectangle const &extent, int width, int height, QgsRasterBlockFeedback *feedback )
809{
810 QgsDebugMsgLevel( QStringLiteral( "extent:\n%1" ).arg( extent.toString() ), 4 );
811 QgsDebugMsgLevel( QStringLiteral( "width = %1 height = %2" ).arg( width ).arg( height ), 4 );
812 if ( !mInput )
813 {
814 QgsDebugMsgLevel( QStringLiteral( "Input not set" ), 4 );
815 return new QgsRasterBlock();
816 }
817
818 if ( feedback && feedback->isCanceled() )
819 return new QgsRasterBlock();
820
821 if ( ! mSrcCRS.isValid() || ! mDestCRS.isValid() || mSrcCRS == mDestCRS )
822 {
823 QgsDebugMsgLevel( QStringLiteral( "No projection necessary" ), 4 );
824 return mInput->block( bandNo, extent, width, height, feedback );
825 }
826
828 const QgsCoordinateTransform inverseCt = mSrcDatumTransform != -1 || mDestDatumTransform != -1 ?
829 QgsCoordinateTransform( mDestCRS, mSrcCRS, mDestDatumTransform, mSrcDatumTransform ) : QgsCoordinateTransform( mDestCRS, mSrcCRS, mTransformContext ) ;
831
832 ProjectorData pd( extent, width, height, mInput, inverseCt, mPrecision, feedback );
833
834 if ( feedback && feedback->isCanceled() )
835 return new QgsRasterBlock();
836
837 QgsDebugMsgLevel( QStringLiteral( "srcExtent:\n%1" ).arg( pd.srcExtent().toString() ), 4 );
838 QgsDebugMsgLevel( QStringLiteral( "srcCols = %1 srcRows = %2" ).arg( pd.srcCols() ).arg( pd.srcRows() ), 4 );
839
840 // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers
841 if ( pd.srcRows() <= 0 || pd.srcCols() <= 0 )
842 {
843 QgsDebugMsgLevel( QStringLiteral( "Zero srcRows or srcCols" ), 4 );
844 return new QgsRasterBlock();
845 }
846
847 std::unique_ptr< QgsRasterBlock > inputBlock( mInput->block( bandNo, pd.srcExtent(), pd.srcCols(), pd.srcRows(), feedback ) );
848 if ( !inputBlock || inputBlock->isEmpty() )
849 {
850 QgsDebugMsg( QStringLiteral( "No raster data!" ) );
851 return new QgsRasterBlock();
852 }
853
854 const qgssize pixelSize = static_cast<qgssize>( QgsRasterBlock::typeSize( mInput->dataType( bandNo ) ) );
855
856 std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock( inputBlock->dataType(), width, height ) );
857 if ( inputBlock->hasNoDataValue() )
858 {
859 outputBlock->setNoDataValue( inputBlock->noDataValue() );
860 }
861 if ( !outputBlock->isValid() )
862 {
863 QgsDebugMsg( QStringLiteral( "Cannot create block" ) );
864 return outputBlock.release();
865 }
866
867 // set output to no data, it should be fast
868 outputBlock->setIsNoData();
869
870 // No data: because isNoData()/setIsNoData() is slow with respect to simple memcpy,
871 // we use if only if necessary:
872 // 1) no data value exists (numerical) -> memcpy, not necessary isNoData()/setIsNoData()
873 // 2) no data value does not exist but it may contain no data (numerical no data bitmap)
874 // -> must use isNoData()/setIsNoData()
875 // 3) no data are not used (no no data value, no no data bitmap) -> simple memcpy
876 // 4) image - simple memcpy
877
878 // To copy no data values stored in bitmaps we have to use isNoData()/setIsNoData(),
879 // we cannot fill output block with no data because we use memcpy for data, not setValue().
880 const bool doNoData = !QgsRasterBlock::typeIsNumeric( inputBlock->dataType() ) && inputBlock->hasNoData() && !inputBlock->hasNoDataValue();
881
882 outputBlock->setIsNoData();
883
884 int srcRow, srcCol;
885 for ( int i = 0; i < height; ++i )
886 {
887 if ( feedback && feedback->isCanceled() )
888 break;
889 for ( int j = 0; j < width; ++j )
890 {
891 const bool inside = pd.srcRowCol( i, j, &srcRow, &srcCol );
892 if ( !inside ) continue; // we have everything set to no data
893
894 const qgssize srcIndex = static_cast< qgssize >( srcRow ) * pd.srcCols() + srcCol;
895
896 // isNoData() may be slow so we check doNoData first
897 if ( doNoData && inputBlock->isNoData( srcRow, srcCol ) )
898 {
899 outputBlock->setIsNoData( i, j );
900 continue;
901 }
902
903 const qgssize destIndex = static_cast< qgssize >( i ) * width + j;
904 char *srcBits = inputBlock->bits( srcIndex );
905 char *destBits = outputBlock->bits( destIndex );
906 if ( !srcBits )
907 {
908 // QgsDebugMsg( QStringLiteral( "Cannot get input block data: row = %1 col = %2" ).arg( i ).arg( j ) );
909 continue;
910 }
911 if ( !destBits )
912 {
913 // QgsDebugMsg( QStringLiteral( "Cannot set output block data: srcRow = %1 srcCol = %2" ).arg( srcRow ).arg( srcCol ) );
914 continue;
915 }
916 memcpy( destBits, srcBits, pixelSize );
917 outputBlock->setIsData( i, j );
918 }
919 }
920
921 return outputBlock.release();
922}
923
924bool QgsRasterProjector::destExtentSize( const QgsRectangle &srcExtent, int srcXSize, int srcYSize,
925 QgsRectangle &destExtent, int &destXSize, int &destYSize )
926{
927 if ( srcExtent.isEmpty() || srcXSize <= 0 || srcYSize <= 0 )
928 {
929 return false;
930 }
931
933 const QgsCoordinateTransform ct = mSrcDatumTransform != -1 || mDestDatumTransform != -1 ?
934 QgsCoordinateTransform( mSrcCRS, mDestCRS, mSrcDatumTransform, mDestDatumTransform ) : QgsCoordinateTransform( mSrcCRS, mDestCRS, mTransformContext ) ;
936
937 return extentSize( ct, srcExtent, srcXSize, srcYSize, destExtent, destXSize, destYSize );
938}
939
941 const QgsRectangle &srcExtent, int srcXSize, int srcYSize,
942 QgsRectangle &destExtent, int &destXSize, int &destYSize )
943{
944 if ( srcExtent.isEmpty() || srcXSize <= 0 || srcYSize <= 0 )
945 {
946 return false;
947 }
948
949 QgsCoordinateTransform extentTransform = ct;
950 extentTransform.setBallparkTransformsAreAppropriate( true );
951
952 destExtent = extentTransform.transformBoundingBox( srcExtent );
953
954 // We reproject pixel rectangle from 9 points matrix of source extent, of course, it gives
955 // bigger xRes,yRes than reprojected edges (envelope)
956 constexpr int steps = 3;
957 const double srcXStep = srcExtent.width() / steps;
958 const double srcYStep = srcExtent.height() / steps;
959 const double srcXRes = srcExtent.width() / srcXSize;
960 const double srcYRes = srcExtent.height() / srcYSize;
961 double destXRes = std::numeric_limits<double>::max();
962 double destYRes = std::numeric_limits<double>::max();
963 double maxXRes = 0;
964 double maxYRes = 0;
965
966 for ( int i = 0; i < steps; i++ )
967 {
968 const double x = srcExtent.xMinimum() + i * srcXStep;
969 for ( int j = 0; j < steps; j++ )
970 {
971 const double y = srcExtent.yMinimum() + j * srcYStep;
972 const QgsRectangle srcRectangle( x - srcXRes / 2, y - srcYRes / 2, x + srcXRes / 2, y + srcYRes / 2 );
973 try
974 {
975 const QgsRectangle destRectangle = extentTransform.transformBoundingBox( srcRectangle );
976 if ( destRectangle.width() > 0 )
977 {
978 destXRes = std::min( destXRes, destRectangle.width() );
979 if ( destRectangle.width() > maxXRes )
980 maxXRes = destRectangle.width();
981 }
982 if ( destRectangle.height() > 0 )
983 {
984 destYRes = std::min( destYRes, destRectangle.height() );
985 if ( destRectangle.height() > maxYRes )
986 maxYRes = destRectangle.height();
987 }
988 }
989 catch ( QgsCsException & )
990 {
991
992 }
993 }
994 }
995
996 // Limit resolution to 1/10th of the maximum resolution to avoid issues
997 // with for example WebMercator at high northings that result in very small
998 // latitude differences.
999 if ( destXRes < 0.1 * maxXRes )
1000 {
1001 destXRes = 0.1 * maxXRes;
1002 }
1003 if ( destYRes < 0.1 * maxYRes )
1004 {
1005 destYRes = 0.1 * maxYRes;
1006 }
1007 if ( destXRes == 0 || destExtent.width() / destXRes > std::numeric_limits<int>::max() )
1008 return false;
1009 if ( destYRes == 0 || destExtent.height() / destYRes > std::numeric_limits<int>::max() )
1010 return false;
1011
1012 destXSize = std::max( 1, static_cast< int >( destExtent.width() / destXRes ) );
1013 destYSize = std::max( 1, static_cast< int >( destExtent.height() / destYRes ) );
1014
1015 return true;
1016}
1017
DataType
Raster data types.
Definition: qgis.h:129
@ UnknownDataType
Unknown or unspecified type.
This class represents a coordinate reference system (CRS).
bool isValid() const
Returns whether this CRS is correctly initialized and usable.
Contains information about the context in which a coordinate transform is executed.
Class for doing transforms between two map coordinate systems.
void setBallparkTransformsAreAppropriate(bool appropriate)
Sets whether approximate "ballpark" results are appropriate for this coordinate transform.
bool isValid() const
Returns true if the coordinate transform is valid, ie both the source and destination CRS have been s...
QgsRectangle transformBoundingBox(const QgsRectangle &rectangle, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool handle180Crossover=false) const SIP_THROW(QgsCsException)
Transforms a rectangle from the source CRS to the destination CRS.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:66
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
A class to represent a 2D point.
Definition: qgspointxy.h:59
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
Definition: qgspointxy.cpp:51
double sqrDist(double x, double y) const SIP_HOLDGIL
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspointxy.h:190
void setX(double x) SIP_HOLDGIL
Sets the x value of the point.
Definition: qgspointxy.h:122
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
void setY(double y) SIP_HOLDGIL
Sets the y value of the point.
Definition: qgspointxy.h:132
Feedback object tailored for raster block reading.
Raster data container.
static bool typeIsNumeric(Qgis::DataType type)
Returns true if data type is numeric.
static int typeSize(Qgis::DataType dataType) SIP_HOLDGIL
Returns the size in bytes for the specified dataType.
Base class for raster data providers.
virtual QgsRasterDataProvider::ProviderCapabilities providerCapabilities() const
Returns flags containing the supported capabilities of the data provider.
QgsRectangle extent() const override=0
Returns the extent of the layer.
@ ProviderHintCanPerformProviderResampling
Provider can perform resampling (to be opposed to post rendering resampling) (since QGIS 3....
Base class for processing filters like renderers, reprojector, resampler etc.
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 Qgis::DataType dataType(int bandNo) const =0
Returns data type for the band specified by number.
@ Size
Original data source size (and thus resolution) is known, it is not always available,...
virtual int xSize() const
Gets raster size.
virtual int bandCount() const =0
Gets number of bands.
QgsRasterInterface * mInput
virtual int capabilities() const
Returns a bitmask containing the supported capabilities.
virtual int ySize() const
virtual QgsRectangle extent() const
Gets the extent of the interface.
virtual QgsRasterInterface * input() const
Current input.
virtual const QgsRasterInterface * sourceInput() const
Gets source / raw input, the first in pipe, usually provider.
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
QgsRasterProjector * clone() const override
Clone itself, create deep copy.
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)
Qgis::DataType dataType(int bandNo) const override
Returns data type for the band specified by number.
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
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.
Precision
Precision defines if each pixel is reprojected or approximate reprojection based on an approximation ...
@ Exact
Exact, precise but slow.
@ Approximate
Approximate (default), fast but possibly inaccurate.
Precision precision() const
int bandCount() const override
Gets number of bands.
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
QString toString(int precision=16) const
Returns a string representation of form xmin,ymin : xmax,ymax Coordinates will be truncated to the sp...
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:469
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
#define Q_NOWARN_DEPRECATED_POP
Definition: qgis.h:3061
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...
Definition: qgis.h:3032
#define Q_NOWARN_DEPRECATED_PUSH
Definition: qgis.h:3060
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
int precision