16
16
***************************************************************************/
17
17
18
18
#include " qgsrasterdataprovider.h"
19
+ #include " qgscrscache.h"
19
20
#include " qgslogger.h"
20
21
#include " qgsrasterprojector.h"
21
22
#include " qgscoordinatetransform.h"
22
23
23
24
QgsRasterProjector::QgsRasterProjector (
24
25
QgsCoordinateReferenceSystem theSrcCRS,
25
26
QgsCoordinateReferenceSystem theDestCRS,
27
+ int theSrcDatumTransform,
28
+ int theDestDatumTransform,
26
29
QgsRectangle theDestExtent,
27
30
int theDestRows, int theDestCols,
28
31
double theMaxSrcXRes, double theMaxSrcYRes,
29
32
QgsRectangle theExtent )
30
33
: QgsRasterInterface( 0 )
31
34
, mSrcCRS( theSrcCRS )
32
35
, mDestCRS( theDestCRS )
33
- , mCoordinateTransform( theDestCRS, theSrcCRS )
36
+ , mSrcDatumTransform( theSrcDatumTransform )
37
+ , mDestDatumTransform( theDestDatumTransform )
34
38
, mDestExtent( theDestExtent )
35
39
, mExtent( theExtent )
36
40
, mDestRows( theDestRows ), mDestCols( theDestCols )
@@ -46,22 +50,46 @@ QgsRasterProjector::QgsRasterProjector(
46
50
QgsRasterProjector::QgsRasterProjector (
47
51
QgsCoordinateReferenceSystem theSrcCRS,
48
52
QgsCoordinateReferenceSystem theDestCRS,
53
+ QgsRectangle theDestExtent,
54
+ int theDestRows, int theDestCols,
49
55
double theMaxSrcXRes, double theMaxSrcYRes,
50
56
QgsRectangle theExtent )
51
57
: QgsRasterInterface( 0 )
52
58
, mSrcCRS( theSrcCRS )
53
59
, mDestCRS( theDestCRS )
54
- , mCoordinateTransform( theDestCRS, theSrcCRS )
60
+ , mSrcDatumTransform( -1 )
61
+ , mDestDatumTransform( -1 )
62
+ , mDestExtent( theDestExtent )
55
63
, mExtent( theExtent )
64
+ , mDestRows( theDestRows ), mDestCols( theDestCols )
56
65
, pHelperTop( 0 ), pHelperBottom( 0 )
57
66
, mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
58
67
{
59
68
QgsDebugMsg ( " Entered" );
69
+ QgsDebugMsg ( " theDestExtent = " + theDestExtent.toString () );
70
+
71
+ calc ();
60
72
}
61
73
62
- QgsRasterProjector::QgsRasterProjector ()
74
+ QgsRasterProjector::QgsRasterProjector (
75
+ QgsCoordinateReferenceSystem theSrcCRS,
76
+ QgsCoordinateReferenceSystem theDestCRS,
77
+ double theMaxSrcXRes, double theMaxSrcYRes,
78
+ QgsRectangle theExtent )
63
79
: QgsRasterInterface( 0 )
80
+ , mSrcCRS( theSrcCRS )
81
+ , mDestCRS( theDestCRS )
82
+ , mSrcDatumTransform( -1 )
83
+ , mDestDatumTransform( -1 )
84
+ , mExtent( theExtent )
64
85
, pHelperTop( 0 ), pHelperBottom( 0 )
86
+ , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
87
+ {
88
+ QgsDebugMsg ( " Entered" );
89
+ }
90
+
91
+ QgsRasterProjector::QgsRasterProjector ()
92
+ : QgsRasterInterface( 0 ), mSrcDatumTransform( -1 ), mDestDatumTransform( -1 ) , pHelperTop( 0 ), pHelperBottom( 0 )
65
93
{
66
94
QgsDebugMsg ( " Entered" );
67
95
}
@@ -71,11 +99,11 @@ QgsRasterProjector::QgsRasterProjector( const QgsRasterProjector &projector )
71
99
{
72
100
mSrcCRS = projector.mSrcCRS ;
73
101
mDestCRS = projector.mDestCRS ;
102
+ mSrcDatumTransform = projector.mSrcDatumTransform ;
103
+ mDestDatumTransform = projector.mDestDatumTransform ;
74
104
mMaxSrcXRes = projector.mMaxSrcXRes ;
75
105
mMaxSrcYRes = projector.mMaxSrcYRes ;
76
106
mExtent = projector.mExtent ;
77
- mCoordinateTransform .setSourceCrs ( mSrcCRS );
78
- mCoordinateTransform .setDestCRS ( mDestCRS );
79
107
}
80
108
81
109
QgsRasterProjector & QgsRasterProjector::operator =( const QgsRasterProjector & projector )
@@ -84,11 +112,11 @@ QgsRasterProjector & QgsRasterProjector::operator=( const QgsRasterProjector & p
84
112
{
85
113
mSrcCRS = projector.mSrcCRS ;
86
114
mDestCRS = projector.mDestCRS ;
115
+ mSrcDatumTransform = projector.mSrcDatumTransform ;
116
+ mDestDatumTransform = projector.mDestDatumTransform ;
87
117
mMaxSrcXRes = projector.mMaxSrcXRes ;
88
118
mMaxSrcYRes = projector.mMaxSrcYRes ;
89
119
mExtent = projector.mExtent ;
90
- mCoordinateTransform .setSourceCrs ( mSrcCRS );
91
- mCoordinateTransform .setDestCRS ( mDestCRS );
92
120
}
93
121
return *this ;
94
122
}
@@ -97,6 +125,8 @@ QgsRasterInterface * QgsRasterProjector::clone() const
97
125
{
98
126
QgsDebugMsg ( " Entered" );
99
127
QgsRasterProjector * projector = new QgsRasterProjector ( mSrcCRS , mDestCRS , mMaxSrcXRes , mMaxSrcYRes , mExtent );
128
+ projector->mSrcDatumTransform = mSrcDatumTransform ;
129
+ projector->mDestDatumTransform = mDestDatumTransform ;
100
130
return projector;
101
131
}
102
132
@@ -120,12 +150,12 @@ QGis::DataType QgsRasterProjector::dataType( int bandNo ) const
120
150
return QGis::UnknownDataType;
121
151
}
122
152
123
- void QgsRasterProjector::setCRS ( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS )
153
+ void QgsRasterProjector::setCRS ( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS, int srcDatumTransform, int destDatumTransform )
124
154
{
125
155
mSrcCRS = theSrcCRS;
126
156
mDestCRS = theDestCRS;
127
- mCoordinateTransform . setSourceCrs ( theDestCRS ) ;
128
- mCoordinateTransform . setDestCRS ( theSrcCRS ) ;
157
+ mSrcDatumTransform = srcDatumTransform ;
158
+ mDestDatumTransform = destDatumTransform ;
129
159
}
130
160
131
161
void QgsRasterProjector::calc ()
@@ -166,6 +196,8 @@ void QgsRasterProjector::calc()
166
196
double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes ;
167
197
mSqrTolerance = myDestRes * myDestRes;
168
198
199
+ const QgsCoordinateTransform* ct = QgsCoordinateTransformCache::instance ()->transform ( mDestCRS .authid (), mSrcCRS .authid (), mDestDatumTransform , mSrcDatumTransform );
200
+
169
201
// Initialize the matrix by corners and middle points
170
202
mCPCols = mCPRows = 3 ;
171
203
for ( int i = 0 ; i < mCPRows ; i++ )
@@ -184,20 +216,20 @@ void QgsRasterProjector::calc()
184
216
}
185
217
for ( int i = 0 ; i < mCPRows ; i++ )
186
218
{
187
- calcRow ( i );
219
+ calcRow ( i, ct );
188
220
}
189
221
190
222
while ( true )
191
223
{
192
- bool myColsOK = checkCols ();
224
+ bool myColsOK = checkCols ( ct );
193
225
if ( !myColsOK )
194
226
{
195
- insertRows ();
227
+ insertRows ( ct );
196
228
}
197
- bool myRowsOK = checkRows ();
229
+ bool myRowsOK = checkRows ( ct );
198
230
if ( !myRowsOK )
199
231
{
200
- insertCols ();
232
+ insertCols ( ct );
201
233
}
202
234
if ( myColsOK && myRowsOK )
203
235
{
@@ -432,19 +464,19 @@ void QgsRasterProjector::nextHelper()
432
464
mHelperTopRow ++;
433
465
}
434
466
435
- void QgsRasterProjector::srcRowCol ( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
467
+ void QgsRasterProjector::srcRowCol ( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
436
468
{
437
469
if ( mApproximate )
438
470
{
439
471
approximateSrcRowCol ( theDestRow, theDestCol, theSrcRow, theSrcCol );
440
472
}
441
473
else
442
474
{
443
- preciseSrcRowCol ( theDestRow, theDestCol, theSrcRow, theSrcCol );
475
+ preciseSrcRowCol ( theDestRow, theDestCol, theSrcRow, theSrcCol, ct );
444
476
}
445
477
}
446
478
447
- void QgsRasterProjector::preciseSrcRowCol ( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
479
+ void QgsRasterProjector::preciseSrcRowCol ( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
448
480
{
449
481
#ifdef QGISDEBUG
450
482
QgsDebugMsgLevel ( QString ( " theDestRow = %1" ).arg ( theDestRow ), 5 );
@@ -460,7 +492,10 @@ void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *
460
492
QgsDebugMsgLevel ( QString ( " x = %1 y = %2" ).arg ( x ).arg ( y ), 5 );
461
493
#endif
462
494
463
- mCoordinateTransform .transformInPlace ( x, y, z );
495
+ if ( ct )
496
+ {
497
+ ct->transformInPlace ( x, y, z );
498
+ }
464
499
465
500
#ifdef QGISDEBUG
466
501
QgsDebugMsgLevel ( QString ( " x = %1 y = %2" ).arg ( x ).arg ( y ), 5 );
@@ -545,7 +580,7 @@ void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, i
545
580
Q_ASSERT ( *theSrcCol < mSrcCols );
546
581
}
547
582
548
- void QgsRasterProjector::insertRows ()
583
+ void QgsRasterProjector::insertRows ( const QgsCoordinateTransform* ct )
549
584
{
550
585
for ( int r = 0 ; r < mCPRows - 1 ; r++ )
551
586
{
@@ -563,11 +598,11 @@ void QgsRasterProjector::insertRows()
563
598
mCPRows += mCPRows - 1 ;
564
599
for ( int r = 1 ; r < mCPRows - 1 ; r += 2 )
565
600
{
566
- calcRow ( r );
601
+ calcRow ( r, ct );
567
602
}
568
603
}
569
604
570
- void QgsRasterProjector::insertCols ()
605
+ void QgsRasterProjector::insertCols ( const QgsCoordinateTransform* ct )
571
606
{
572
607
for ( int r = 0 ; r < mCPRows ; r++ )
573
608
{
@@ -582,20 +617,27 @@ void QgsRasterProjector::insertCols()
582
617
mCPCols += mCPCols - 1 ;
583
618
for ( int c = 1 ; c < mCPCols - 1 ; c += 2 )
584
619
{
585
- calcCol ( c );
620
+ calcCol ( c, ct );
586
621
}
587
622
588
623
}
589
624
590
- void QgsRasterProjector::calcCP ( int theRow, int theCol )
625
+ void QgsRasterProjector::calcCP ( int theRow, int theCol, const QgsCoordinateTransform* ct )
591
626
{
592
627
double myDestX, myDestY;
593
628
destPointOnCPMatrix ( theRow, theCol, &myDestX, &myDestY );
594
629
QgsPoint myDestPoint ( myDestX, myDestY );
595
630
try
596
631
{
597
- mCPMatrix [theRow][theCol] = mCoordinateTransform .transform ( myDestPoint );
598
- mCPLegalMatrix [theRow][theCol] = true ;
632
+ if ( ct )
633
+ {
634
+ mCPMatrix [theRow][theCol] = ct->transform ( myDestPoint );
635
+ mCPLegalMatrix [theRow][theCol] = true ;
636
+ }
637
+ else
638
+ {
639
+ mCPLegalMatrix [theRow][theCol] = false ;
640
+ }
599
641
}
600
642
catch ( QgsCsException &e )
601
643
{
@@ -605,30 +647,35 @@ void QgsRasterProjector::calcCP( int theRow, int theCol )
605
647
}
606
648
}
607
649
608
- bool QgsRasterProjector::calcRow ( int theRow )
650
+ bool QgsRasterProjector::calcRow ( int theRow, const QgsCoordinateTransform* ct )
609
651
{
610
652
QgsDebugMsgLevel ( QString ( " theRow = %1" ).arg ( theRow ), 3 );
611
653
for ( int i = 0 ; i < mCPCols ; i++ )
612
654
{
613
- calcCP ( theRow, i );
655
+ calcCP ( theRow, i, ct );
614
656
}
615
657
616
658
return true ;
617
659
}
618
660
619
- bool QgsRasterProjector::calcCol ( int theCol )
661
+ bool QgsRasterProjector::calcCol ( int theCol, const QgsCoordinateTransform* ct )
620
662
{
621
663
QgsDebugMsgLevel ( QString ( " theCol = %1" ).arg ( theCol ), 3 );
622
664
for ( int i = 0 ; i < mCPRows ; i++ )
623
665
{
624
- calcCP ( i, theCol );
666
+ calcCP ( i, theCol, ct );
625
667
}
626
668
627
669
return true ;
628
670
}
629
671
630
- bool QgsRasterProjector::checkCols ()
672
+ bool QgsRasterProjector::checkCols ( const QgsCoordinateTransform* ct )
631
673
{
674
+ if ( !ct )
675
+ {
676
+ return false ;
677
+ }
678
+
632
679
for ( int c = 0 ; c < mCPCols ; c++ )
633
680
{
634
681
for ( int r = 1 ; r < mCPRows - 1 ; r += 2 )
@@ -649,7 +696,7 @@ bool QgsRasterProjector::checkCols()
649
696
}
650
697
try
651
698
{
652
- QgsPoint myDestApprox = mCoordinateTransform . transform ( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
699
+ QgsPoint myDestApprox = ct-> transform ( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
653
700
double mySqrDist = myDestApprox.sqrDist ( myDestPoint );
654
701
if ( mySqrDist > mSqrTolerance )
655
702
{
@@ -667,8 +714,13 @@ bool QgsRasterProjector::checkCols()
667
714
return true ;
668
715
}
669
716
670
- bool QgsRasterProjector::checkRows ()
717
+ bool QgsRasterProjector::checkRows ( const QgsCoordinateTransform* ct )
671
718
{
719
+ if ( !ct )
720
+ {
721
+ return false ;
722
+ }
723
+
672
724
for ( int r = 0 ; r < mCPRows ; r++ )
673
725
{
674
726
for ( int c = 1 ; c < mCPCols - 1 ; c += 2 )
@@ -689,7 +741,7 @@ bool QgsRasterProjector::checkRows()
689
741
}
690
742
try
691
743
{
692
- QgsPoint myDestApprox = mCoordinateTransform . transform ( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
744
+ QgsPoint myDestApprox = ct-> transform ( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
693
745
double mySqrDist = myDestApprox.sqrDist ( myDestPoint );
694
746
if ( mySqrDist > mSqrTolerance )
695
747
{
@@ -775,12 +827,18 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex
775
827
// we cannot fill output block with no data because we use memcpy for data, not setValue().
776
828
bool doNoData = inputBlock->hasNoData () && !inputBlock->hasNoDataValue ();
777
829
830
+ const QgsCoordinateTransform* ct = 0 ;
831
+ if ( !mApproximate )
832
+ {
833
+ ct = QgsCoordinateTransformCache::instance ()->transform ( mDestCRS .authid (), mSrcCRS .authid (), mDestDatumTransform , mSrcDatumTransform );
834
+ }
835
+
778
836
int srcRow, srcCol;
779
837
for ( int i = 0 ; i < height; ++i )
780
838
{
781
839
for ( int j = 0 ; j < width; ++j )
782
840
{
783
- srcRowCol ( i, j, &srcRow, &srcCol );
841
+ srcRowCol ( i, j, &srcRow, &srcCol, ct );
784
842
qgssize srcIndex = ( qgssize )srcRow * mSrcCols + srcCol;
785
843
QgsDebugMsgLevel ( QString ( " row = %1 col = %2 srcRow = %3 srcCol = %4" ).arg ( i ).arg ( j ).arg ( srcRow ).arg ( srcCol ), 5 );
786
844
0 commit comments