Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalignraster.cpp
3 : : --------------------------------------
4 : : Date : June 2015
5 : : Copyright : (C) 2015 by Martin Dobias
6 : : Email : wonder dot sk at gmail dot com
7 : : ***************************************************************************
8 : : * *
9 : : * This program is free software; you can redistribute it and/or modify *
10 : : * it under the terms of the GNU General Public License as published by *
11 : : * the Free Software Foundation; either version 2 of the License, or *
12 : : * (at your option) any later version. *
13 : : * *
14 : : ***************************************************************************/
15 : :
16 : : #include "qgsalignraster.h"
17 : :
18 : : #include <gdalwarper.h>
19 : : #include <ogr_srs_api.h>
20 : : #include <cpl_conv.h>
21 : : #include <limits>
22 : :
23 : : #include <QPair>
24 : : #include <QString>
25 : :
26 : : #include "qgscoordinatereferencesystem.h"
27 : : #include "qgsrectangle.h"
28 : :
29 : :
30 : 0 : static double ceil_with_tolerance( double value )
31 : : {
32 : 0 : if ( std::fabs( value - std::round( value ) ) < 1e-6 )
33 : 0 : return std::round( value );
34 : : else
35 : 0 : return std::ceil( value );
36 : 0 : }
37 : :
38 : 0 : static double floor_with_tolerance( double value )
39 : : {
40 : 0 : if ( std::fabs( value - std::round( value ) ) < 1e-6 )
41 : 0 : return std::round( value );
42 : : else
43 : 0 : return std::floor( value );
44 : 0 : }
45 : :
46 : 0 : static double fmod_with_tolerance( double num, double denom )
47 : : {
48 : 0 : return num - floor_with_tolerance( num / denom ) * denom;
49 : : }
50 : :
51 : :
52 : 0 : static QgsRectangle transform_to_extent( const double *geotransform, double xSize, double ySize )
53 : : {
54 : 0 : QgsRectangle r( geotransform[0],
55 : 0 : geotransform[3],
56 : 0 : geotransform[0] + geotransform[1] * xSize,
57 : 0 : geotransform[3] + geotransform[5] * ySize );
58 : 0 : r.normalize();
59 : 0 : return r;
60 : : }
61 : :
62 : :
63 : 0 : static int CPL_STDCALL _progress( double dfComplete, const char *pszMessage, void *pProgressArg )
64 : : {
65 : : Q_UNUSED( pszMessage )
66 : :
67 : 0 : QgsAlignRaster::ProgressHandler *handler = ( ( QgsAlignRaster * ) pProgressArg )->progressHandler();
68 : 0 : if ( handler )
69 : 0 : return handler->progress( dfComplete );
70 : : else
71 : 0 : return true;
72 : 0 : }
73 : :
74 : :
75 : 0 : static CPLErr rescalePreWarpChunkProcessor( void *pKern, void *pArg )
76 : : {
77 : 0 : GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
78 : 0 : double cellsize = ( ( double * )pArg )[0];
79 : :
80 : 0 : for ( int nBand = 0; nBand < kern->nBands; ++nBand )
81 : : {
82 : 0 : float *bandData = ( float * ) kern->papabySrcImage[nBand];
83 : 0 : for ( int nLine = 0; nLine < kern->nSrcYSize; ++nLine )
84 : : {
85 : 0 : for ( int nPixel = 0; nPixel < kern->nSrcXSize; ++nPixel )
86 : : {
87 : 0 : bandData[nPixel] /= cellsize;
88 : 0 : }
89 : 0 : bandData += kern->nSrcXSize;
90 : 0 : }
91 : 0 : }
92 : 0 : return CE_None;
93 : : }
94 : :
95 : 0 :
96 : 0 : static CPLErr rescalePostWarpChunkProcessor( void *pKern, void *pArg )
97 : 0 : {
98 : 0 : GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
99 : 0 : double cellsize = ( ( double * )pArg )[1];
100 : :
101 : 0 : for ( int nBand = 0; nBand < kern->nBands; ++nBand )
102 : : {
103 : 0 : float *bandData = ( float * ) kern->papabyDstImage[nBand];
104 : 0 : for ( int nLine = 0; nLine < kern->nDstYSize; ++nLine )
105 : : {
106 : 0 : for ( int nPixel = 0; nPixel < kern->nDstXSize; ++nPixel )
107 : : {
108 : 0 : bandData[nPixel] *= cellsize;
109 : 0 : }
110 : 0 : bandData += kern->nDstXSize;
111 : 0 : }
112 : 0 : }
113 : 0 : return CE_None;
114 : : }
115 : :
116 : :
117 : :
118 : 0 : QgsAlignRaster::QgsAlignRaster()
119 : : {
120 : : // parameters
121 : 0 : mCellSizeX = mCellSizeY = 0;
122 : 0 : mGridOffsetX = mGridOffsetY = 0;
123 : 0 : mClipExtent[0] = mClipExtent[1] = mClipExtent[2] = mClipExtent[3] = 0;
124 : :
125 : : // derived variables
126 : 0 : mXSize = mYSize = 0;
127 : 0 : for ( int i = 0; i < 6; ++i )
128 : 0 : mGeoTransform[i] = 0;
129 : 0 : }
130 : :
131 : 0 : void QgsAlignRaster::setClipExtent( double xmin, double ymin, double xmax, double ymax )
132 : : {
133 : 0 : mClipExtent[0] = xmin;
134 : 0 : mClipExtent[1] = ymin;
135 : 0 : mClipExtent[2] = xmax;
136 : 0 : mClipExtent[3] = ymax;
137 : 0 : }
138 : :
139 : 0 : void QgsAlignRaster::setClipExtent( const QgsRectangle &extent )
140 : : {
141 : 0 : setClipExtent( extent.xMinimum(), extent.yMinimum(),
142 : 0 : extent.xMaximum(), extent.yMaximum() );
143 : 0 : }
144 : :
145 : 0 : QgsRectangle QgsAlignRaster::clipExtent() const
146 : : {
147 : 0 : return QgsRectangle( mClipExtent[0], mClipExtent[1],
148 : 0 : mClipExtent[2], mClipExtent[3] );
149 : : }
150 : :
151 : :
152 : 0 : bool QgsAlignRaster::setParametersFromRaster( const QString &filename, const QString &destWkt, QSizeF customCellSize, QPointF customGridOffset )
153 : : {
154 : 0 : return setParametersFromRaster( RasterInfo( filename ), destWkt, customCellSize, customGridOffset );
155 : 0 : }
156 : :
157 : 0 : bool QgsAlignRaster::setParametersFromRaster( const RasterInfo &rasterInfo, const QString &customCRSWkt, QSizeF customCellSize, QPointF customGridOffset )
158 : : {
159 : 0 : if ( customCRSWkt.isEmpty() || customCRSWkt == rasterInfo.crs() )
160 : : {
161 : : // use ref. layer to init input
162 : 0 : mCrsWkt = rasterInfo.crs();
163 : :
164 : 0 : if ( !customCellSize.isValid() )
165 : : {
166 : 0 : mCellSizeX = rasterInfo.cellSize().width();
167 : 0 : mCellSizeY = rasterInfo.cellSize().height();
168 : 0 : }
169 : : else
170 : : {
171 : 0 : mCellSizeX = customCellSize.width();
172 : 0 : mCellSizeY = customCellSize.height();
173 : : }
174 : :
175 : 0 : if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
176 : : {
177 : 0 : if ( !customCellSize.isValid() )
178 : : {
179 : : // using original raster's grid offset to be aligned with origin
180 : 0 : mGridOffsetX = rasterInfo.gridOffset().x();
181 : 0 : mGridOffsetY = rasterInfo.gridOffset().y();
182 : 0 : }
183 : : else
184 : : {
185 : : // if using custom cell size: offset so that we are aligned
186 : : // with the original raster's origin point
187 : 0 : mGridOffsetX = fmod_with_tolerance( rasterInfo.origin().x(), customCellSize.width() );
188 : 0 : mGridOffsetY = fmod_with_tolerance( rasterInfo.origin().y(), customCellSize.height() );
189 : : }
190 : 0 : }
191 : : else
192 : : {
193 : 0 : mGridOffsetX = customGridOffset.x();
194 : 0 : mGridOffsetY = customGridOffset.y();
195 : : }
196 : 0 : }
197 : : else
198 : : {
199 : 0 : QSizeF cs;
200 : 0 : QPointF go;
201 : 0 : if ( !suggestedWarpOutput( rasterInfo, customCRSWkt, &cs, &go ) )
202 : : {
203 : 0 : mCrsWkt = QStringLiteral( "_error_" );
204 : 0 : mCellSizeX = mCellSizeY = 0;
205 : 0 : mGridOffsetX = mGridOffsetY = 0;
206 : 0 : return false;
207 : : }
208 : :
209 : 0 : mCrsWkt = customCRSWkt;
210 : :
211 : 0 : if ( !customCellSize.isValid() )
212 : : {
213 : 0 : mCellSizeX = cs.width();
214 : 0 : mCellSizeY = cs.height();
215 : 0 : }
216 : : else
217 : : {
218 : 0 : mCellSizeX = customCellSize.width();
219 : 0 : mCellSizeY = customCellSize.height();
220 : : }
221 : :
222 : 0 : if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
223 : : {
224 : 0 : mGridOffsetX = go.x();
225 : 0 : mGridOffsetY = go.y();
226 : 0 : }
227 : : else
228 : : {
229 : 0 : mGridOffsetX = customGridOffset.x();
230 : 0 : mGridOffsetY = customGridOffset.x();
231 : : }
232 : : }
233 : 0 : return true;
234 : 0 : }
235 : :
236 : :
237 : 0 : bool QgsAlignRaster::checkInputParameters()
238 : : {
239 : 0 : mErrorMessage.clear();
240 : :
241 : 0 : if ( mCrsWkt == QLatin1String( "_error_" ) )
242 : : {
243 : 0 : mErrorMessage = QObject::tr( "Unable to reproject." );
244 : 0 : return false;
245 : : }
246 : :
247 : 0 : if ( mCellSizeX == 0 || mCellSizeY == 0 )
248 : : {
249 : 0 : mErrorMessage = QObject::tr( "Cell size must not be zero." );
250 : 0 : return false;
251 : : }
252 : :
253 : 0 : mXSize = mYSize = 0;
254 : 0 : std::fill_n( mGeoTransform, 6, 0 );
255 : :
256 : 0 : double finalExtent[4] = { 0, 0, 0, 0 };
257 : :
258 : : // for each raster: determine their extent in projected cfg
259 : 0 : for ( int i = 0; i < mRasters.count(); ++i )
260 : : {
261 : 0 : Item &r = mRasters[i];
262 : :
263 : 0 : RasterInfo info( r.inputFilename );
264 : :
265 : 0 : QSizeF cs;
266 : 0 : QgsRectangle extent;
267 : 0 : if ( !suggestedWarpOutput( info, mCrsWkt, &cs, nullptr, &extent ) )
268 : : {
269 : 0 : mErrorMessage = QString( "Failed to get suggested warp output.\n\n"
270 : : "File:\n%1\n\n"
271 : : "Source WKT:\n%2\n\nDestination WKT:\n%3" )
272 : 0 : .arg( r.inputFilename,
273 : 0 : info.mCrsWkt,
274 : 0 : mCrsWkt );
275 : 0 : return false;
276 : : }
277 : :
278 : 0 : r.srcCellSizeInDestCRS = cs.width() * cs.height();
279 : :
280 : 0 : if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
281 : : {
282 : 0 : // initialize with the first layer
283 : 0 : finalExtent[0] = extent.xMinimum();
284 : 0 : finalExtent[1] = extent.yMinimum();
285 : 0 : finalExtent[2] = extent.xMaximum();
286 : 0 : finalExtent[3] = extent.yMaximum();
287 : 0 : }
288 : : else
289 : : {
290 : : // use intersection of rects
291 : 0 : if ( extent.xMinimum() > finalExtent[0] ) finalExtent[0] = extent.xMinimum();
292 : 0 : if ( extent.yMinimum() > finalExtent[1] ) finalExtent[1] = extent.yMinimum();
293 : 0 : if ( extent.xMaximum() < finalExtent[2] ) finalExtent[2] = extent.xMaximum();
294 : 0 : if ( extent.yMaximum() < finalExtent[3] ) finalExtent[3] = extent.yMaximum();
295 : : }
296 : 0 : }
297 : :
298 : : // count in extra clip extent (if present)
299 : : // 1. align requested rect to the grid - extend the rect if necessary
300 : : // 2. intersect with clip extent with final extent
301 : :
302 : 0 : if ( !( mClipExtent[0] == 0 && mClipExtent[1] == 0 && mClipExtent[2] == 0 && mClipExtent[3] == 0 ) )
303 : : {
304 : : // extend clip extent to grid
305 : 0 : double clipX0 = floor_with_tolerance( ( mClipExtent[0] - mGridOffsetX ) / mCellSizeX ) * mCellSizeX + mGridOffsetX;
306 : 0 : double clipY0 = floor_with_tolerance( ( mClipExtent[1] - mGridOffsetY ) / mCellSizeY ) * mCellSizeY + mGridOffsetY;
307 : 0 : double clipX1 = ceil_with_tolerance( ( mClipExtent[2] - clipX0 ) / mCellSizeX ) * mCellSizeX + clipX0;
308 : 0 : double clipY1 = ceil_with_tolerance( ( mClipExtent[3] - clipY0 ) / mCellSizeY ) * mCellSizeY + clipY0;
309 : 0 : if ( clipX0 > finalExtent[0] ) finalExtent[0] = clipX0;
310 : 0 : if ( clipY0 > finalExtent[1] ) finalExtent[1] = clipY0;
311 : 0 : if ( clipX1 < finalExtent[2] ) finalExtent[2] = clipX1;
312 : 0 : if ( clipY1 < finalExtent[3] ) finalExtent[3] = clipY1;
313 : 0 : }
314 : :
315 : : // align to grid - shrink the rect if necessary
316 : : // output raster grid configuration (with no rotation/shear)
317 : : // ... and raster width/height
318 : :
319 : 0 : double originX = ceil_with_tolerance( ( finalExtent[0] - mGridOffsetX ) / mCellSizeX ) * mCellSizeX + mGridOffsetX;
320 : 0 : double originY = ceil_with_tolerance( ( finalExtent[1] - mGridOffsetY ) / mCellSizeY ) * mCellSizeY + mGridOffsetY;
321 : 0 : int xSize = floor_with_tolerance( ( finalExtent[2] - originX ) / mCellSizeX );
322 : 0 : int ySize = floor_with_tolerance( ( finalExtent[3] - originY ) / mCellSizeY );
323 : :
324 : 0 : if ( xSize <= 0 || ySize <= 0 )
325 : : {
326 : 0 : mErrorMessage = QObject::tr( "No common intersecting area." );
327 : 0 : return false;
328 : : }
329 : :
330 : 0 : mXSize = xSize;
331 : 0 : mYSize = ySize;
332 : :
333 : : // build final geotransform...
334 : 0 : mGeoTransform[0] = originX;
335 : 0 : mGeoTransform[1] = mCellSizeX;
336 : 0 : mGeoTransform[2] = 0;
337 : 0 : mGeoTransform[3] = originY + ( mCellSizeY * ySize );
338 : 0 : mGeoTransform[4] = 0;
339 : 0 : mGeoTransform[5] = -mCellSizeY;
340 : :
341 : 0 : return true;
342 : 0 : }
343 : :
344 : :
345 : 0 : QSize QgsAlignRaster::alignedRasterSize() const
346 : : {
347 : 0 : return QSize( mXSize, mYSize );
348 : : }
349 : :
350 : 0 : QgsRectangle QgsAlignRaster::alignedRasterExtent() const
351 : : {
352 : 0 : return transform_to_extent( mGeoTransform, mXSize, mYSize );
353 : : }
354 : :
355 : :
356 : 0 : bool QgsAlignRaster::run()
357 : : {
358 : 0 : mErrorMessage.clear();
359 : :
360 : : // consider extent of all layers and setup geotransform and output grid size
361 : 0 : if ( !checkInputParameters() )
362 : 0 : return false;
363 : :
364 : : //dump();
365 : :
366 : 0 : const auto constMRasters = mRasters;
367 : 0 : for ( const Item &r : constMRasters )
368 : : {
369 : 0 : if ( !createAndWarp( r ) )
370 : 0 : return false;
371 : : }
372 : 0 : return true;
373 : 0 : }
374 : :
375 : :
376 : 0 : void QgsAlignRaster::dump() const
377 : : {
378 : 0 : qDebug( "---ALIGN------------------" );
379 : 0 : qDebug( "wkt %s", mCrsWkt.toLatin1().constData() );
380 : 0 : qDebug( "w/h %d,%d", mXSize, mYSize );
381 : 0 : qDebug( "transform" );
382 : 0 : qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[0], mGeoTransform[1], mGeoTransform[2] );
383 : 0 : qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[3], mGeoTransform[4], mGeoTransform[5] );
384 : :
385 : 0 : QgsRectangle e = transform_to_extent( mGeoTransform, mXSize, mYSize );
386 : 0 : qDebug( "extent %s", e.toString().toLatin1().constData() );
387 : 0 : }
388 : :
389 : 0 : int QgsAlignRaster::suggestedReferenceLayer() const
390 : : {
391 : 0 : int bestIndex = -1;
392 : 0 : double bestCellArea = INFINITY;
393 : 0 : QSizeF cs;
394 : 0 : int i = 0;
395 : :
396 : : // using WGS84 as a destination CRS... but maybe some projected CRS
397 : : // would be a better a choice to more accurately compute areas?
398 : : // (Why earth is not flat???)
399 : 0 : QgsCoordinateReferenceSystem destCRS( QStringLiteral( "EPSG:4326" ) );
400 : 0 : QString destWkt = destCRS.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL );
401 : :
402 : 0 : const auto constMRasters = mRasters;
403 : 0 : for ( const Item &raster : constMRasters )
404 : : {
405 : 0 : if ( !suggestedWarpOutput( RasterInfo( raster.inputFilename ), destWkt, &cs ) )
406 : 0 : return false;
407 : :
408 : 0 : double cellArea = cs.width() * cs.height();
409 : 0 : if ( cellArea < bestCellArea )
410 : : {
411 : 0 : bestCellArea = cellArea;
412 : 0 : bestIndex = i;
413 : 0 : }
414 : 0 : ++i;
415 : : }
416 : :
417 : 0 : return bestIndex;
418 : 0 : }
419 : :
420 : :
421 : 0 : bool QgsAlignRaster::createAndWarp( const Item &raster )
422 : : {
423 : 0 : GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );
424 : 0 : if ( !hDriver )
425 : : {
426 : 0 : mErrorMessage = QStringLiteral( "GDALGetDriverByName(GTiff) failed." );
427 : 0 : return false;
428 : : }
429 : :
430 : : // Open the source file.
431 : 0 : gdal::dataset_unique_ptr hSrcDS( GDALOpen( raster.inputFilename.toLocal8Bit().constData(), GA_ReadOnly ) );
432 : 0 : if ( !hSrcDS )
433 : : {
434 : 0 : mErrorMessage = QObject::tr( "Unable to open input file: %1" ).arg( raster.inputFilename );
435 : 0 : return false;
436 : : }
437 : :
438 : : // Create output with same datatype as first input band.
439 : :
440 : 0 : int bandCount = GDALGetRasterCount( hSrcDS.get() );
441 : 0 : GDALDataType eDT = GDALGetRasterDataType( GDALGetRasterBand( hSrcDS.get(), 1 ) );
442 : :
443 : : // Create the output file.
444 : 0 : gdal::dataset_unique_ptr hDstDS( GDALCreate( hDriver, raster.outputFilename.toLocal8Bit().constData(), mXSize, mYSize,
445 : 0 : bandCount, eDT, nullptr ) );
446 : 0 : if ( !hDstDS )
447 : : {
448 : 0 : mErrorMessage = QObject::tr( "Unable to create output file: %1" ).arg( raster.outputFilename );
449 : 0 : return false;
450 : : }
451 : :
452 : : // Write out the projection definition.
453 : 0 : GDALSetProjection( hDstDS.get(), mCrsWkt.toLatin1().constData() );
454 : 0 : GDALSetGeoTransform( hDstDS.get(), mGeoTransform );
455 : :
456 : : // Copy the color table, if required.
457 : 0 : GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand( hSrcDS.get(), 1 ) );
458 : 0 : if ( hCT )
459 : 0 : GDALSetRasterColorTable( GDALGetRasterBand( hDstDS.get(), 1 ), hCT );
460 : :
461 : : // -----------------------------------------------------------------------
462 : :
463 : : // Setup warp options.
464 : 0 : gdal::warp_options_unique_ptr psWarpOptions( GDALCreateWarpOptions() );
465 : 0 : psWarpOptions->hSrcDS = hSrcDS.get();
466 : 0 : psWarpOptions->hDstDS = hDstDS.get();
467 : :
468 : 0 : psWarpOptions->nBandCount = GDALGetRasterCount( hSrcDS.get() );
469 : 0 : psWarpOptions->panSrcBands = ( int * ) CPLMalloc( sizeof( int ) * psWarpOptions->nBandCount );
470 : 0 : psWarpOptions->panDstBands = ( int * ) CPLMalloc( sizeof( int ) * psWarpOptions->nBandCount );
471 : 0 : for ( int i = 0; i < psWarpOptions->nBandCount; ++i )
472 : : {
473 : 0 : psWarpOptions->panSrcBands[i] = i + 1;
474 : 0 : psWarpOptions->panDstBands[i] = i + 1;
475 : 0 : }
476 : :
477 : 0 : psWarpOptions->eResampleAlg = static_cast< GDALResampleAlg >( raster.resampleMethod );
478 : :
479 : : // our progress function
480 : 0 : psWarpOptions->pfnProgress = _progress;
481 : 0 : psWarpOptions->pProgressArg = this;
482 : :
483 : : // Establish reprojection transformer.
484 : 0 : psWarpOptions->pTransformerArg =
485 : 0 : GDALCreateGenImgProjTransformer( hSrcDS.get(), GDALGetProjectionRef( hSrcDS.get() ),
486 : 0 : hDstDS.get(), GDALGetProjectionRef( hDstDS.get() ),
487 : : FALSE, 0.0, 1 );
488 : 0 : psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
489 : :
490 : : double rescaleArg[2];
491 : 0 : if ( raster.rescaleValues )
492 : : {
493 : 0 : rescaleArg[0] = raster.srcCellSizeInDestCRS; // source cell size
494 : 0 : rescaleArg[1] = mCellSizeX * mCellSizeY; // destination cell size
495 : 0 : psWarpOptions->pfnPreWarpChunkProcessor = rescalePreWarpChunkProcessor;
496 : 0 : psWarpOptions->pfnPostWarpChunkProcessor = rescalePostWarpChunkProcessor;
497 : 0 : psWarpOptions->pPreWarpProcessorArg = rescaleArg;
498 : 0 : psWarpOptions->pPostWarpProcessorArg = rescaleArg;
499 : : // force use of float32 data type as that is what our pre/post-processor uses
500 : 0 : psWarpOptions->eWorkingDataType = GDT_Float32;
501 : 0 : }
502 : :
503 : : // Initialize and execute the warp operation.
504 : 0 : GDALWarpOperation oOperation;
505 : 0 : oOperation.Initialize( psWarpOptions.get() );
506 : 0 : oOperation.ChunkAndWarpImage( 0, 0, mXSize, mYSize );
507 : :
508 : 0 : GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
509 : 0 : return true;
510 : 0 : }
511 : :
512 : 0 : bool QgsAlignRaster::suggestedWarpOutput( const QgsAlignRaster::RasterInfo &info, const QString &destWkt, QSizeF *cellSize, QPointF *gridOffset, QgsRectangle *rect )
513 : : {
514 : : // Create a transformer that maps from source pixel/line coordinates
515 : : // to destination georeferenced coordinates (not destination
516 : : // pixel line). We do that by omitting the destination dataset
517 : : // handle (setting it to nullptr).
518 : 0 : void *hTransformArg = GDALCreateGenImgProjTransformer( info.mDataset.get(), info.mCrsWkt.toLatin1().constData(), nullptr, destWkt.toLatin1().constData(), FALSE, 0, 1 );
519 : 0 : if ( !hTransformArg )
520 : 0 : return false;
521 : :
522 : : // Get approximate output georeferenced bounds and resolution for file.
523 : : double adfDstGeoTransform[6];
524 : : double extents[4];
525 : 0 : int nPixels = 0, nLines = 0;
526 : : CPLErr eErr;
527 : 0 : eErr = GDALSuggestedWarpOutput2( info.mDataset.get(),
528 : 0 : GDALGenImgProjTransform, hTransformArg,
529 : 0 : adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
530 : 0 : GDALDestroyGenImgProjTransformer( hTransformArg );
531 : :
532 : 0 : if ( eErr != CE_None )
533 : 0 : return false;
534 : :
535 : 0 : QSizeF cs( std::fabs( adfDstGeoTransform[1] ), std::fabs( adfDstGeoTransform[5] ) );
536 : :
537 : 0 : if ( rect )
538 : 0 : *rect = QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
539 : 0 : if ( cellSize )
540 : 0 : *cellSize = cs;
541 : 0 : if ( gridOffset )
542 : 0 : *gridOffset = QPointF( fmod_with_tolerance( adfDstGeoTransform[0], cs.width() ),
543 : 0 : fmod_with_tolerance( adfDstGeoTransform[3], cs.height() ) );
544 : 0 : return true;
545 : 0 : }
546 : :
547 : :
548 : : //----------
549 : :
550 : :
551 : 0 : QgsAlignRaster::RasterInfo::RasterInfo( const QString &layerpath )
552 : : {
553 : 0 : mDataset.reset( GDALOpen( layerpath.toLocal8Bit().constData(), GA_ReadOnly ) );
554 : 0 : if ( !mDataset )
555 : 0 : return;
556 : :
557 : 0 : mXSize = GDALGetRasterXSize( mDataset.get() );
558 : 0 : mYSize = GDALGetRasterYSize( mDataset.get() );
559 : :
560 : 0 : ( void ) GDALGetGeoTransform( mDataset.get(), mGeoTransform );
561 : :
562 : : // TODO: may be null or empty string
563 : 0 : mCrsWkt = QString::fromLatin1( GDALGetProjectionRef( mDataset.get() ) );
564 : :
565 : 0 : mBandCnt = GDALGetBandNumber( mDataset.get() );
566 : 0 : }
567 : :
568 : 0 : QSizeF QgsAlignRaster::RasterInfo::cellSize() const
569 : : {
570 : 0 : return QSizeF( std::fabs( mGeoTransform[1] ), std::fabs( mGeoTransform[5] ) );
571 : : }
572 : :
573 : 0 : QPointF QgsAlignRaster::RasterInfo::gridOffset() const
574 : : {
575 : 0 : return QPointF( fmod_with_tolerance( mGeoTransform[0], cellSize().width() ),
576 : 0 : fmod_with_tolerance( mGeoTransform[3], cellSize().height() ) );
577 : : }
578 : :
579 : 0 : QgsRectangle QgsAlignRaster::RasterInfo::extent() const
580 : : {
581 : 0 : return transform_to_extent( mGeoTransform, mXSize, mYSize );
582 : : }
583 : :
584 : 0 : QPointF QgsAlignRaster::RasterInfo::origin() const
585 : : {
586 : 0 : return QPointF( mGeoTransform[0], mGeoTransform[3] );
587 : : }
588 : :
589 : 0 : void QgsAlignRaster::RasterInfo::dump() const
590 : : {
591 : 0 : qDebug( "---RASTER INFO------------------" );
592 : 0 : qDebug( "wkt %s", mCrsWkt.toLatin1().constData() );
593 : 0 : qDebug( "w/h %d,%d", mXSize, mYSize );
594 : 0 : qDebug( "cell x/y %f,%f", cellSize().width(), cellSize().width() );
595 : :
596 : 0 : QgsRectangle r = extent();
597 : 0 : qDebug( "extent %s", r.toString().toLatin1().constData() );
598 : :
599 : 0 : qDebug( "transform" );
600 : 0 : qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[0], mGeoTransform[1], mGeoTransform[2] );
601 : 0 : qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[3], mGeoTransform[4], mGeoTransform[5] );
602 : 0 : }
603 : :
604 : 0 : double QgsAlignRaster::RasterInfo::identify( double mx, double my )
605 : : {
606 : 0 : GDALRasterBandH hBand = GDALGetRasterBand( mDataset.get(), 1 );
607 : :
608 : : // must not be rotated in order for this to work
609 : 0 : int px = int( ( mx - mGeoTransform[0] ) / mGeoTransform[1] );
610 : 0 : int py = int( ( my - mGeoTransform[3] ) / mGeoTransform[5] );
611 : :
612 : 0 : float *pafScanline = ( float * ) CPLMalloc( sizeof( float ) );
613 : 0 : CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,
614 : 0 : pafScanline, 1, 1, GDT_Float32, 0, 0 );
615 : 0 : double value = err == CE_None ? pafScanline[0] : std::numeric_limits<double>::quiet_NaN();
616 : 0 : CPLFree( pafScanline );
617 : :
618 : 0 : return value;
619 : : }
620 : :
|
