1 /***************************************************************************
2 qgsgeoreftransform.cpp - Encapsulates GCP-based parameter estimation and
3 reprojection for different transformation models.
4 --------------------------------------
5 Date : 18-Feb-2009
6 Copyright : (c) 2009 by Manuel Massing
7 Email : m.massing at warped-space.de
8 ***************************************************************************
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 ***************************************************************************/
16
17 #include "qgsgeoreftransform.h"
18
19 #include <gdal.h>
20 #include <gdal_alg.h>
21
22 #include <cmath>
23
24 #include <cassert>
25 #include <limits>
26
QgsGeorefTransform(const QgsGeorefTransform & other)27 QgsGeorefTransform::QgsGeorefTransform( const QgsGeorefTransform &other )
28 {
29 selectTransformParametrisation( other.mTransformParametrisation );
30 }
31
QgsGeorefTransform(TransformMethod parametrisation)32 QgsGeorefTransform::QgsGeorefTransform( TransformMethod parametrisation )
33 {
34 selectTransformParametrisation( parametrisation );
35 }
36
37 QgsGeorefTransform::QgsGeorefTransform() = default;
38
39 QgsGeorefTransform::~QgsGeorefTransform() = default;
40
transformParametrisation() const41 QgsGeorefTransform::TransformMethod QgsGeorefTransform::transformParametrisation() const
42 {
43 return mTransformParametrisation;
44 }
45
selectTransformParametrisation(TransformMethod parametrisation)46 void QgsGeorefTransform::selectTransformParametrisation( TransformMethod parametrisation )
47 {
48 if ( parametrisation != mTransformParametrisation )
49 {
50 mGeorefTransformImplementation.reset( QgsGcpTransformerInterface::create( parametrisation ) );
51 mParametersInitialized = false;
52 mTransformParametrisation = parametrisation;
53 }
54 }
55
setRasterChangeCoords(const QString & fileRaster)56 void QgsGeorefTransform::setRasterChangeCoords( const QString &fileRaster )
57 {
58 mRasterChangeCoords.setRaster( fileRaster );
59 }
60
providesAccurateInverseTransformation() const61 bool QgsGeorefTransform::providesAccurateInverseTransformation() const
62 {
63 return ( mTransformParametrisation == TransformMethod::Linear
64 || mTransformParametrisation == TransformMethod::Helmert
65 || mTransformParametrisation == TransformMethod::PolynomialOrder1 );
66 }
67
parametersInitialized() const68 bool QgsGeorefTransform::parametersInitialized() const
69 {
70 return mParametersInitialized;
71 }
72
clone() const73 QgsGcpTransformerInterface *QgsGeorefTransform::clone() const
74 {
75 std::unique_ptr< QgsGeorefTransform > res( new QgsGeorefTransform( *this ) );
76 res->updateParametersFromGcps( mSourceCoordinates, mDestinationCoordinates, mInvertYAxis );
77 return res.release();
78 }
79
updateParametersFromGcps(const QVector<QgsPointXY> & sourceCoordinates,const QVector<QgsPointXY> & destinationCoordinates,bool invertYAxis)80 bool QgsGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
81 {
82 mSourceCoordinates = sourceCoordinates;
83 mDestinationCoordinates = destinationCoordinates;
84 mInvertYAxis = invertYAxis;
85
86 if ( !mGeorefTransformImplementation )
87 {
88 return false;
89 }
90 if ( sourceCoordinates.size() != destinationCoordinates.size() ) // Defensive sanity check
91 {
92 throw ( std::domain_error( "Internal error: GCP mapping is not one-to-one" ) );
93 }
94 if ( sourceCoordinates.size() < minimumGcpCount() )
95 {
96 return false;
97 }
98 if ( mRasterChangeCoords.hasCrs() )
99 {
100 const QVector<QgsPointXY> pixelCoordsCorrect = mRasterChangeCoords.getPixelCoords( sourceCoordinates );
101 mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( sourceCoordinates, pixelCoordsCorrect, invertYAxis );
102 }
103 else
104 {
105 mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( sourceCoordinates, destinationCoordinates, invertYAxis );
106 }
107 return mParametersInitialized;
108 }
109
minimumGcpCount() const110 int QgsGeorefTransform::minimumGcpCount() const
111 {
112 return mGeorefTransformImplementation ? mGeorefTransformImplementation->minimumGcpCount() : 0;
113 }
114
method() const115 QgsGcpTransformerInterface::TransformMethod QgsGeorefTransform::method() const
116 {
117 return mGeorefTransformImplementation ? mGeorefTransformImplementation->method() : TransformMethod::InvalidTransform;
118 }
119
GDALTransformer() const120 GDALTransformerFunc QgsGeorefTransform::GDALTransformer() const
121 {
122 return mGeorefTransformImplementation ? mGeorefTransformImplementation->GDALTransformer() : nullptr;
123 }
124
GDALTransformerArgs() const125 void *QgsGeorefTransform::GDALTransformerArgs() const
126 {
127 return mGeorefTransformImplementation ? mGeorefTransformImplementation->GDALTransformerArgs() : nullptr;
128 }
129
transformRasterToWorld(const QgsPointXY & raster,QgsPointXY & world)130 bool QgsGeorefTransform::transformRasterToWorld( const QgsPointXY &raster, QgsPointXY &world )
131 {
132 // flip y coordinate due to different CS orientation
133 const QgsPointXY raster_flipped( raster.x(), -raster.y() );
134 return gdal_transform( raster_flipped, world, 0 );
135 }
136
transformWorldToRaster(const QgsPointXY & world,QgsPointXY & raster)137 bool QgsGeorefTransform::transformWorldToRaster( const QgsPointXY &world, QgsPointXY &raster )
138 {
139 const bool success = gdal_transform( world, raster, 1 );
140 // flip y coordinate due to different CS orientation
141 raster.setY( -raster.y() );
142 return success;
143 }
144
transform(const QgsPointXY & src,QgsPointXY & dst,bool rasterToWorld)145 bool QgsGeorefTransform::transform( const QgsPointXY &src, QgsPointXY &dst, bool rasterToWorld )
146 {
147 return rasterToWorld ? transformRasterToWorld( src, dst ) : transformWorldToRaster( src, dst );
148 }
149
getLinearOriginScale(QgsPointXY & origin,double & scaleX,double & scaleY) const150 bool QgsGeorefTransform::getLinearOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const
151 {
152 if ( transformParametrisation() != TransformMethod::Linear )
153 {
154 return false;
155 }
156 if ( !mGeorefTransformImplementation || !parametersInitialized() )
157 {
158 return false;
159 }
160 QgsLinearGeorefTransform *transform = dynamic_cast<QgsLinearGeorefTransform *>( mGeorefTransformImplementation.get() );
161 return transform && transform->getOriginScale( origin, scaleX, scaleY );
162 }
163
getOriginScaleRotation(QgsPointXY & origin,double & scaleX,double & scaleY,double & rotation) const164 bool QgsGeorefTransform::getOriginScaleRotation( QgsPointXY &origin, double &scaleX, double &scaleY, double &rotation ) const
165 {
166
167 if ( mTransformParametrisation == TransformMethod::Linear )
168 {
169 rotation = 0.0;
170 QgsLinearGeorefTransform *transform = dynamic_cast<QgsLinearGeorefTransform *>( mGeorefTransformImplementation.get() );
171 return transform && transform->getOriginScale( origin, scaleX, scaleY );
172 }
173 else if ( mTransformParametrisation == TransformMethod::Helmert )
174 {
175 double scale;
176 QgsHelmertGeorefTransform *transform = dynamic_cast<QgsHelmertGeorefTransform *>( mGeorefTransformImplementation.get() );
177 if ( !transform || ! transform->getOriginScaleRotation( origin, scale, rotation ) )
178 {
179 return false;
180 }
181 scaleX = scale;
182 scaleY = scale;
183 return true;
184 }
185 return false;
186 }
187
188
gdal_transform(const QgsPointXY & src,QgsPointXY & dst,int dstToSrc) const189 bool QgsGeorefTransform::gdal_transform( const QgsPointXY &src, QgsPointXY &dst, int dstToSrc ) const
190 {
191 // Copy the source coordinate for inplace transform
192 double x = src.x();
193 double y = src.y();
194
195 if ( !QgsGcpTransformerInterface::transform( x, y, dstToSrc == 1 ) )
196 return false;
197
198 dst.setX( x );
199 dst.setY( y );
200 return true;
201 }
202
203
204