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