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