1 : /******************************************************************************
2 : * $Id: vrtwarped.cpp 17852 2009-10-18 11:15:09Z rouault $
3 : *
4 : * Project: Virtual GDAL Datasets
5 : * Purpose: Implementation of VRTWarpedRasterBand *and VRTWarpedDataset.
6 : * Author: Frank Warmerdam <warmerdam@pobox.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "vrtdataset.h"
31 : #include "cpl_minixml.h"
32 : #include "cpl_string.h"
33 : #include "gdalwarper.h"
34 : #include <cassert>
35 :
36 : CPL_CVSID("$Id: vrtwarped.cpp 17852 2009-10-18 11:15:09Z rouault $");
37 :
38 : /************************************************************************/
39 : /* GDALAutoCreateWarpedVRT() */
40 : /************************************************************************/
41 :
42 : /**
43 : * Create virtual warped dataset automatically.
44 : *
45 : * This function will create a warped virtual file representing the
46 : * input image warped into the target coordinate system. A GenImgProj
47 : * transformation is created to accomplish any required GCP/Geotransform
48 : * warp and reprojection to the target coordinate system. The output virtual
49 : * dataset will be "northup" in the target coordinate system. The
50 : * GDALSuggestedWarpOutput() function is used to determine the bounds and
51 : * resolution of the output virtual file which should be large enough to
52 : * include all the input image
53 : *
54 : * Note that the constructed GDALDatasetH will acquire one or more references
55 : * to the passed in hSrcDS. Reference counting semantics on the source
56 : * dataset should be honoured. That is, don't just GDALClose() it unless it
57 : * was opened with GDALOpenShared().
58 : *
59 : * The returned dataset will have no associated filename for itself. If you
60 : * want to write the virtual dataset description to a file, use the
61 : * GDALSetDescription() function (or SetDescription() method) on the dataset
62 : * to assign a filename before it is closed.
63 : *
64 : * @param hSrcDS The source dataset.
65 : *
66 : * @param pszSrcWKT The coordinate system of the source image. If NULL, it
67 : * will be read from the source image.
68 : *
69 : * @param pszDstWKT The coordinate system to convert to. If NULL no change
70 : * of coordinate system will take place.
71 : *
72 : * @param eResampleAlg One of GRA_NearestNeighbour, GRA_Bilinear, GRA_Cubic or
73 : * GRA_CubicSpline. Controls the sampling method used.
74 : *
75 : * @param dfMaxError Maximum error measured in input pixels that is allowed in
76 : * approximating the transformation (0.0 for exact calculations).
77 : *
78 : * @param psOptions Additional warp options, normally NULL.
79 : *
80 : * @return NULL on failure, or a new virtual dataset handle on success.
81 : */
82 :
83 : GDALDatasetH CPL_STDCALL
84 1 : GDALAutoCreateWarpedVRT( GDALDatasetH hSrcDS,
85 : const char *pszSrcWKT,
86 : const char *pszDstWKT,
87 : GDALResampleAlg eResampleAlg,
88 : double dfMaxError,
89 : const GDALWarpOptions *psOptionsIn )
90 :
91 : {
92 : GDALWarpOptions *psWO;
93 : int i;
94 :
95 1 : VALIDATE_POINTER1( hSrcDS, "GDALAutoCreateWarpedVRT", NULL );
96 :
97 : /* -------------------------------------------------------------------- */
98 : /* Populate the warp options. */
99 : /* -------------------------------------------------------------------- */
100 1 : if( psOptionsIn != NULL )
101 0 : psWO = GDALCloneWarpOptions( psOptionsIn );
102 : else
103 1 : psWO = GDALCreateWarpOptions();
104 :
105 1 : psWO->eResampleAlg = eResampleAlg;
106 :
107 1 : psWO->hSrcDS = hSrcDS;
108 :
109 1 : psWO->nBandCount = GDALGetRasterCount( hSrcDS );
110 1 : psWO->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
111 1 : psWO->panDstBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
112 :
113 4 : for( i = 0; i < psWO->nBandCount; i++ )
114 : {
115 3 : psWO->panSrcBands[i] = i+1;
116 3 : psWO->panDstBands[i] = i+1;
117 : }
118 :
119 : /* TODO: should fill in no data where available */
120 :
121 : /* -------------------------------------------------------------------- */
122 : /* Create the transformer. */
123 : /* -------------------------------------------------------------------- */
124 1 : psWO->pfnTransformer = GDALGenImgProjTransform;
125 : psWO->pTransformerArg =
126 : GDALCreateGenImgProjTransformer( psWO->hSrcDS, pszSrcWKT,
127 : NULL, pszDstWKT,
128 1 : TRUE, 1.0, 0 );
129 :
130 1 : if( psWO->pTransformerArg == NULL )
131 : {
132 0 : GDALDestroyWarpOptions( psWO );
133 0 : return NULL;
134 : }
135 :
136 : /* -------------------------------------------------------------------- */
137 : /* Figure out the desired output bounds and resolution. */
138 : /* -------------------------------------------------------------------- */
139 : double adfDstGeoTransform[6];
140 : int nDstPixels, nDstLines;
141 : CPLErr eErr;
142 :
143 : eErr =
144 : GDALSuggestedWarpOutput( hSrcDS, psWO->pfnTransformer,
145 : psWO->pTransformerArg,
146 1 : adfDstGeoTransform, &nDstPixels, &nDstLines );
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Update the transformer to include an output geotransform */
150 : /* back to pixel/line coordinates. */
151 : /* */
152 : /* -------------------------------------------------------------------- */
153 : GDALSetGenImgProjTransformerDstGeoTransform(
154 1 : psWO->pTransformerArg, adfDstGeoTransform );
155 :
156 : /* -------------------------------------------------------------------- */
157 : /* Do we want to apply an approximating transformation? */
158 : /* -------------------------------------------------------------------- */
159 1 : if( dfMaxError > 0.0 )
160 : {
161 : psWO->pTransformerArg =
162 : GDALCreateApproxTransformer( psWO->pfnTransformer,
163 : psWO->pTransformerArg,
164 0 : dfMaxError );
165 0 : psWO->pfnTransformer = GDALApproxTransform;
166 : }
167 :
168 : /* -------------------------------------------------------------------- */
169 : /* Create the VRT file. */
170 : /* -------------------------------------------------------------------- */
171 : GDALDatasetH hDstDS;
172 :
173 : hDstDS = GDALCreateWarpedVRT( hSrcDS, nDstPixels, nDstLines,
174 1 : adfDstGeoTransform, psWO );
175 :
176 1 : GDALDestroyWarpOptions( psWO );
177 :
178 1 : if( pszDstWKT != NULL )
179 0 : GDALSetProjection( hDstDS, pszDstWKT );
180 1 : else if( pszSrcWKT != NULL )
181 0 : GDALSetProjection( hDstDS, pszDstWKT );
182 1 : else if( GDALGetGCPCount( hSrcDS ) > 0 )
183 1 : GDALSetProjection( hDstDS, GDALGetGCPProjection( hSrcDS ) );
184 : else
185 0 : GDALSetProjection( hDstDS, GDALGetProjectionRef( hSrcDS ) );
186 :
187 1 : return hDstDS;
188 : }
189 :
190 : /************************************************************************/
191 : /* GDALCreateWarpedVRT() */
192 : /************************************************************************/
193 :
194 : /**
195 : * Create virtual warped dataset.
196 : *
197 : * This function will create a warped virtual file representing the
198 : * input image warped based on a provided transformation. Output bounds
199 : * and resolution are provided explicitly.
200 : *
201 : * Note that the constructed GDALDatasetH will acquire one or more references
202 : * to the passed in hSrcDS. Reference counting semantics on the source
203 : * dataset should be honoured. That is, don't just GDALClose() it unless it
204 : * was opened with GDALOpenShared().
205 : *
206 : * @param hSrcDS The source dataset.
207 : *
208 : * @param nPixels Width of the virtual warped dataset to create
209 : *
210 : * @param nLines Height of the virtual warped dataset to create
211 : *
212 : * @param padfGeoTransform Geotransform matrix of the virtual warped dataset to create
213 : *
214 : * @param psOptions Warp options. Must be different from NULL.
215 : *
216 : * @return NULL on failure, or a new virtual dataset handle on success.
217 : */
218 :
219 : GDALDatasetH CPL_STDCALL
220 1 : GDALCreateWarpedVRT( GDALDatasetH hSrcDS,
221 : int nPixels, int nLines, double *padfGeoTransform,
222 : GDALWarpOptions *psOptions )
223 :
224 : {
225 1 : VALIDATE_POINTER1( hSrcDS, "GDALCreateWarpedVRT", NULL );
226 :
227 : /* -------------------------------------------------------------------- */
228 : /* Create the VRTDataset and populate it with bands. */
229 : /* -------------------------------------------------------------------- */
230 1 : VRTWarpedDataset *poDS = new VRTWarpedDataset( nPixels, nLines );
231 : int i;
232 :
233 1 : psOptions->hDstDS = (GDALDatasetH) poDS;
234 :
235 1 : poDS->SetGeoTransform( padfGeoTransform );
236 :
237 4 : for( i = 0; i < psOptions->nBandCount; i++ )
238 : {
239 : VRTWarpedRasterBand *poBand;
240 : GDALRasterBand *poSrcBand = (GDALRasterBand *)
241 3 : GDALGetRasterBand( hSrcDS, i+1 );
242 :
243 3 : poDS->AddBand( poSrcBand->GetRasterDataType(), NULL );
244 :
245 3 : poBand = (VRTWarpedRasterBand *) poDS->GetRasterBand( i+1 );
246 3 : poBand->CopyCommonInfoFrom( poSrcBand );
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Initialize the warp on the VRTWarpedDataset. */
251 : /* -------------------------------------------------------------------- */
252 1 : poDS->Initialize( psOptions );
253 :
254 1 : return (GDALDatasetH) poDS;
255 : }
256 :
257 : /************************************************************************/
258 : /* ==================================================================== */
259 : /* VRTWarpedDataset */
260 : /* ==================================================================== */
261 : /************************************************************************/
262 :
263 : /************************************************************************/
264 : /* VRTWarpedDataset() */
265 : /************************************************************************/
266 :
267 39 : VRTWarpedDataset::VRTWarpedDataset( int nXSize, int nYSize )
268 39 : : VRTDataset( nXSize, nYSize )
269 :
270 : {
271 39 : poWarper = NULL;
272 39 : nBlockXSize = 512;
273 39 : nBlockYSize = 128;
274 39 : eAccess = GA_Update;
275 :
276 39 : nOverviewCount = 0;
277 39 : papoOverviews = NULL;
278 39 : }
279 :
280 : /************************************************************************/
281 : /* ~VRTWarpedDataset() */
282 : /************************************************************************/
283 :
284 78 : VRTWarpedDataset::~VRTWarpedDataset()
285 :
286 : {
287 39 : FlushCache();
288 :
289 : /* -------------------------------------------------------------------- */
290 : /* Cleanup overviews. */
291 : /* -------------------------------------------------------------------- */
292 : int iOverview;
293 :
294 41 : for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
295 : {
296 2 : GDALDatasetH hDS = (GDALDatasetH) papoOverviews[iOverview];
297 :
298 2 : if( GDALDereferenceDataset( hDS ) < 1 )
299 : {
300 2 : GDALReferenceDataset( hDS );
301 2 : GDALClose( hDS );
302 : }
303 : }
304 :
305 39 : CPLFree( papoOverviews );
306 :
307 : /* -------------------------------------------------------------------- */
308 : /* Cleanup warper if one is in effect. */
309 : /* -------------------------------------------------------------------- */
310 39 : if( poWarper != NULL )
311 : {
312 39 : const GDALWarpOptions *psWO = poWarper->GetOptions();
313 :
314 : /* -------------------------------------------------------------------- */
315 : /* We take care to only call GDALClose() on psWO->hSrcDS if the */
316 : /* reference count drops to zero. This is makes it so that we */
317 : /* can operate reference counting semantics more-or-less */
318 : /* properly even if the dataset isn't open in shared mode, */
319 : /* though we require that the caller also honour the reference */
320 : /* counting semantics even though it isn't a shared dataset. */
321 : /* -------------------------------------------------------------------- */
322 39 : if( psWO->hSrcDS != NULL )
323 : {
324 39 : if( GDALDereferenceDataset( psWO->hSrcDS ) < 1 )
325 : {
326 35 : GDALReferenceDataset( psWO->hSrcDS );
327 35 : GDALClose( psWO->hSrcDS );
328 : }
329 : }
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* We are responsible for cleaning up the transformer outselves. */
333 : /* -------------------------------------------------------------------- */
334 39 : if( psWO->pTransformerArg != NULL )
335 39 : GDALDestroyTransformer( psWO->pTransformerArg );
336 :
337 39 : delete poWarper;
338 : }
339 78 : }
340 :
341 : /************************************************************************/
342 : /* Initialize() */
343 : /* */
344 : /* Initialize a dataset from passed in warp options. */
345 : /************************************************************************/
346 :
347 5 : CPLErr VRTWarpedDataset::Initialize( void *psWO )
348 :
349 : {
350 5 : if( poWarper != NULL )
351 0 : delete poWarper;
352 :
353 5 : poWarper = new GDALWarpOperation();
354 :
355 : // The act of initializing this warped dataset with this warp options
356 : // will result in our assuming ownership of a reference to the
357 : // hSrcDS.
358 :
359 5 : if( ((GDALWarpOptions *) psWO)->hSrcDS != NULL )
360 5 : GDALReferenceDataset( ((GDALWarpOptions *) psWO)->hSrcDS );
361 :
362 5 : return poWarper->Initialize( (GDALWarpOptions *) psWO );
363 : }
364 :
365 : /************************************************************************/
366 : /* GetFileList() */
367 : /************************************************************************/
368 :
369 0 : char** VRTWarpedDataset::GetFileList()
370 : {
371 0 : char** papszFileList = GDALDataset::GetFileList();
372 :
373 0 : if( poWarper != NULL )
374 : {
375 0 : const GDALWarpOptions *psWO = poWarper->GetOptions();
376 : const char* pszFilename;
377 :
378 0 : if( psWO->hSrcDS != NULL &&
379 : (pszFilename =
380 0 : ((GDALDataset*)psWO->hSrcDS)->GetDescription()) != NULL )
381 : {
382 : VSIStatBufL sStat;
383 0 : if( VSIStatL( pszFilename, &sStat ) == 0 )
384 : {
385 0 : papszFileList = CSLAddString(papszFileList, pszFilename);
386 : }
387 : }
388 : }
389 :
390 0 : return papszFileList;
391 : }
392 :
393 : /************************************************************************/
394 : /* VRTWarpedOverviewTransform() */
395 : /************************************************************************/
396 :
397 : typedef struct {
398 : GDALTransformerInfo sTI;
399 :
400 : GDALTransformerFunc pfnBaseTransformer;
401 : void *pBaseTransformerArg;
402 :
403 : double dfXOverviewFactor;
404 : double dfYOverviewFactor;
405 : } VWOTInfo;
406 :
407 : static
408 258 : int VRTWarpedOverviewTransform( void *pTransformArg, int bDstToSrc,
409 : int nPointCount,
410 : double *padfX, double *padfY, double *padfZ,
411 : int *panSuccess )
412 :
413 : {
414 258 : VWOTInfo *psInfo = (VWOTInfo *) pTransformArg;
415 : int i, bSuccess;
416 :
417 258 : if( bDstToSrc )
418 : {
419 131498 : for( i = 0; i < nPointCount; i++ )
420 : {
421 131240 : padfX[i] *= psInfo->dfXOverviewFactor;
422 131240 : padfY[i] *= psInfo->dfYOverviewFactor;
423 : }
424 : }
425 :
426 : bSuccess = psInfo->pfnBaseTransformer( psInfo->pBaseTransformerArg,
427 : bDstToSrc,
428 : nPointCount, padfX, padfY, padfZ,
429 258 : panSuccess );
430 :
431 258 : if( !bDstToSrc )
432 : {
433 0 : for( i = 0; i < nPointCount; i++ )
434 : {
435 0 : padfX[i] /= psInfo->dfXOverviewFactor;
436 0 : padfY[i] /= psInfo->dfYOverviewFactor;
437 : }
438 : }
439 :
440 258 : return bSuccess;
441 : }
442 :
443 2 : static void VRTWarpedOverviewCleanup(void* pTransformArg)
444 : {
445 2 : VWOTInfo *psInfo = (VWOTInfo *) pTransformArg;
446 :
447 2 : CPLFree( psInfo );
448 2 : }
449 :
450 : /************************************************************************/
451 : /* BuildOverviews() */
452 : /* */
453 : /* For overviews, we actually just build a whole new dataset */
454 : /* with an extra layer of transformation on the warper used to */
455 : /* accomplish downsampling by the desired factor. */
456 : /************************************************************************/
457 :
458 : CPLErr
459 2 : VRTWarpedDataset::IBuildOverviews( const char *pszResampling,
460 : int nOverviews, int *panOverviewList,
461 : int nListBands, int *panBandList,
462 : GDALProgressFunc pfnProgress,
463 : void * pProgressData )
464 :
465 : {
466 : /* -------------------------------------------------------------------- */
467 : /* Initial progress result. */
468 : /* -------------------------------------------------------------------- */
469 2 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
470 : {
471 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
472 0 : return CE_Failure;
473 : }
474 :
475 : /* -------------------------------------------------------------------- */
476 : /* Establish which of the overview levels we already have, and */
477 : /* which are new. */
478 : /* -------------------------------------------------------------------- */
479 2 : int i, nNewOverviews, *panNewOverviewList = NULL;
480 :
481 2 : nNewOverviews = 0;
482 2 : panNewOverviewList = (int *) CPLCalloc(sizeof(int),nOverviews);
483 4 : for( i = 0; i < nOverviews; i++ )
484 : {
485 : int j;
486 :
487 2 : for( j = 0; j < nOverviewCount; j++ )
488 : {
489 : int nOvFactor;
490 0 : VRTWarpedDataset *poOverview = papoOverviews[j];
491 :
492 : nOvFactor = (int)
493 0 : (0.5+GetRasterXSize() / (double) poOverview->GetRasterXSize());
494 :
495 0 : if( nOvFactor == panOverviewList[i]
496 : || nOvFactor == GDALOvLevelAdjust( panOverviewList[i],
497 : GetRasterXSize() ) )
498 0 : panOverviewList[i] *= -1;
499 : }
500 :
501 2 : if( panOverviewList[i] > 0 )
502 2 : panNewOverviewList[nNewOverviews++] = panOverviewList[i];
503 : }
504 :
505 : /* -------------------------------------------------------------------- */
506 : /* Create each missing overview (we don't need to do anything */
507 : /* to update existing overviews). */
508 : /* -------------------------------------------------------------------- */
509 4 : for( i = 0; i < nNewOverviews; i++ )
510 : {
511 : int nOXSize, nOYSize, iBand;
512 : VWOTInfo *psInfo;
513 : VRTWarpedDataset *poOverviewDS;
514 :
515 : /* -------------------------------------------------------------------- */
516 : /* What size should this overview be. */
517 : /* -------------------------------------------------------------------- */
518 2 : nOXSize = (GetRasterXSize() + panNewOverviewList[i] - 1)
519 4 : / panNewOverviewList[i];
520 :
521 2 : nOYSize = (GetRasterYSize() + panNewOverviewList[i] - 1)
522 4 : / panNewOverviewList[i];
523 :
524 : /* -------------------------------------------------------------------- */
525 : /* Create the overview dataset. */
526 : /* -------------------------------------------------------------------- */
527 2 : poOverviewDS = new VRTWarpedDataset( nOXSize, nOYSize );
528 :
529 8 : for( iBand = 0; iBand < GetRasterCount(); iBand++ )
530 : {
531 2 : GDALRasterBand *poOldBand = GetRasterBand(iBand+1);
532 : VRTWarpedRasterBand *poNewBand =
533 : new VRTWarpedRasterBand( poOverviewDS, iBand+1,
534 2 : poOldBand->GetRasterDataType() );
535 :
536 2 : poNewBand->CopyCommonInfoFrom( poOldBand );
537 2 : poOverviewDS->SetBand( iBand+1, poNewBand );
538 : }
539 :
540 2 : nOverviewCount++;
541 : papoOverviews = (VRTWarpedDataset **)
542 2 : CPLRealloc( papoOverviews, sizeof(void*) * nOverviewCount );
543 :
544 2 : papoOverviews[nOverviewCount-1] = poOverviewDS;
545 :
546 : /* -------------------------------------------------------------------- */
547 : /* Prepare update transformation information that will apply */
548 : /* the overview decimation. */
549 : /* -------------------------------------------------------------------- */
550 2 : GDALWarpOptions *psWO = (GDALWarpOptions *) poWarper->GetOptions();
551 2 : psInfo = (VWOTInfo *) CPLCalloc(sizeof(VWOTInfo),1);
552 :
553 2 : strcpy( psInfo->sTI.szSignature, "GTI" );
554 2 : psInfo->sTI.pszClassName = "VRTWarpedOverviewTransform";
555 2 : psInfo->sTI.pfnTransform = VRTWarpedOverviewTransform;
556 2 : psInfo->sTI.pfnCleanup = VRTWarpedOverviewCleanup;
557 2 : psInfo->sTI.pfnSerialize = NULL;
558 :
559 2 : psInfo->pfnBaseTransformer = psWO->pfnTransformer;
560 2 : psInfo->pBaseTransformerArg = psWO->pTransformerArg;
561 :
562 2 : psInfo->dfXOverviewFactor = GetRasterXSize() / (double) nOXSize;
563 2 : psInfo->dfYOverviewFactor = GetRasterYSize() / (double) nOYSize;
564 :
565 : /* -------------------------------------------------------------------- */
566 : /* Initialize the new dataset with adjusted warp options, and */
567 : /* then restore to original condition. */
568 : /* -------------------------------------------------------------------- */
569 2 : psWO->pfnTransformer = VRTWarpedOverviewTransform;
570 2 : psWO->pTransformerArg = psInfo;
571 :
572 2 : poOverviewDS->Initialize( psWO );
573 :
574 2 : psWO->pfnTransformer = psInfo->pfnBaseTransformer;
575 2 : psWO->pTransformerArg = psInfo->pBaseTransformerArg;
576 : }
577 :
578 2 : CPLFree( panNewOverviewList );
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* Progress finished. */
582 : /* -------------------------------------------------------------------- */
583 2 : pfnProgress( 1.0, NULL, pProgressData );
584 :
585 2 : SetNeedsFlush();
586 :
587 2 : return CE_None;
588 : }
589 :
590 : /************************************************************************/
591 : /* GDALInitializeWarpedVRT() */
592 : /************************************************************************/
593 :
594 : /**
595 : * Set warp info on virtual warped dataset.
596 : *
597 : * Initializes all the warping information for a virtual warped dataset.
598 : *
599 : * This method is the same as the C++ method VRTWarpedDataset::Initialize().
600 : *
601 : * @param hDS dataset previously created with the VRT driver, and a
602 : * SUBCLASS of "VRTWarpedDataset".
603 : *
604 : * @param psWO the warp options to apply. Note that ownership of the
605 : * transformation information is taken over by the function though everything
606 : * else remains the property of the caller.
607 : *
608 : * @return CE_None on success or CE_Failure if an error occurs.
609 : */
610 :
611 : CPLErr CPL_STDCALL
612 2 : GDALInitializeWarpedVRT( GDALDatasetH hDS, GDALWarpOptions *psWO )
613 :
614 : {
615 2 : VALIDATE_POINTER1( hDS, "GDALInitializeWarpedVRT", CE_Failure );
616 :
617 2 : return ((VRTWarpedDataset *) hDS)->Initialize( psWO );
618 : }
619 :
620 : /************************************************************************/
621 : /* XMLInit() */
622 : /************************************************************************/
623 :
624 34 : CPLErr VRTWarpedDataset::XMLInit( CPLXMLNode *psTree, const char *pszVRTPath )
625 :
626 : {
627 : CPLErr eErr;
628 :
629 : /* -------------------------------------------------------------------- */
630 : /* Initialize blocksize before calling sub-init so that the */
631 : /* band initializers can get it from the dataset object when */
632 : /* they are created. */
633 : /* -------------------------------------------------------------------- */
634 34 : nBlockXSize = atoi(CPLGetXMLValue(psTree,"BlockXSize","512"));
635 34 : nBlockYSize = atoi(CPLGetXMLValue(psTree,"BlockYSize","128"));
636 :
637 : /* -------------------------------------------------------------------- */
638 : /* Initialize all the general VRT stuff. This will even */
639 : /* create the VRTWarpedRasterBands and initialize them. */
640 : /* -------------------------------------------------------------------- */
641 34 : eErr = VRTDataset::XMLInit( psTree, pszVRTPath );
642 :
643 34 : if( eErr != CE_None )
644 0 : return eErr;
645 :
646 : /* -------------------------------------------------------------------- */
647 : /* Find the GDALWarpOptions XML tree. */
648 : /* -------------------------------------------------------------------- */
649 : CPLXMLNode *psOptionsTree;
650 34 : psOptionsTree = CPLGetXMLNode( psTree, "GDALWarpOptions" );
651 34 : if( psOptionsTree == NULL )
652 : {
653 : CPLError( CE_Failure, CPLE_AppDefined,
654 0 : "Count not find required GDALWarpOptions in XML." );
655 0 : return CE_Failure;
656 : }
657 :
658 : /* -------------------------------------------------------------------- */
659 : /* Adjust the SourceDataset in the warp options to take into */
660 : /* account that it is relative to the VRT if appropriate. */
661 : /* -------------------------------------------------------------------- */
662 : int bRelativeToVRT =
663 : atoi(CPLGetXMLValue(psOptionsTree,
664 34 : "SourceDataset.relativeToVRT", "0" ));
665 :
666 : const char *pszRelativePath = CPLGetXMLValue(psOptionsTree,
667 34 : "SourceDataset", "" );
668 : char *pszAbsolutePath;
669 :
670 34 : if( bRelativeToVRT )
671 : pszAbsolutePath =
672 : CPLStrdup(CPLProjectRelativeFilename( pszVRTPath,
673 32 : pszRelativePath ) );
674 : else
675 2 : pszAbsolutePath = CPLStrdup(pszRelativePath);
676 :
677 34 : CPLSetXMLValue( psOptionsTree, "SourceDataset", pszAbsolutePath );
678 34 : CPLFree( pszAbsolutePath );
679 :
680 : /* -------------------------------------------------------------------- */
681 : /* And instantiate the warp options, and corresponding warp */
682 : /* operation. */
683 : /* -------------------------------------------------------------------- */
684 : GDALWarpOptions *psWO;
685 :
686 34 : psWO = GDALDeserializeWarpOptions( psOptionsTree );
687 34 : if( psWO == NULL )
688 0 : return CE_Failure;
689 :
690 34 : this->eAccess = GA_Update;
691 34 : psWO->hDstDS = this;
692 :
693 : /* -------------------------------------------------------------------- */
694 : /* Instantiate the warp operation. */
695 : /* -------------------------------------------------------------------- */
696 34 : poWarper = new GDALWarpOperation();
697 :
698 34 : eErr = poWarper->Initialize( psWO );
699 34 : if( eErr != CE_None)
700 : {
701 : /* -------------------------------------------------------------------- */
702 : /* We are responsible for cleaning up the transformer outselves. */
703 : /* -------------------------------------------------------------------- */
704 0 : if( psWO->pTransformerArg != NULL )
705 0 : GDALDestroyTransformer( psWO->pTransformerArg );
706 : }
707 :
708 34 : GDALDestroyWarpOptions( psWO );
709 34 : if( eErr != CE_None )
710 : {
711 0 : delete poWarper;
712 0 : poWarper = NULL;
713 : }
714 :
715 : /* -------------------------------------------------------------------- */
716 : /* Generate overviews, if appropriate. */
717 : /* -------------------------------------------------------------------- */
718 : char **papszTokens = CSLTokenizeString(
719 34 : CPLGetXMLValue( psTree, "OverviewList", "" ));
720 : int iOverview;
721 :
722 70 : for( iOverview = 0;
723 35 : papszTokens != NULL && papszTokens[iOverview] != NULL;
724 : iOverview++ )
725 : {
726 1 : int nOvFactor = atoi(papszTokens[iOverview]);
727 :
728 1 : if (nOvFactor > 0)
729 1 : BuildOverviews( "NEAREST", 1, &nOvFactor, 0, NULL, NULL, NULL );
730 : else
731 : CPLError(CE_Failure, CPLE_AppDefined,
732 0 : "Bad value for overview factor : %s", papszTokens[iOverview]);
733 : }
734 :
735 34 : CSLDestroy( papszTokens );
736 :
737 34 : return eErr;
738 : }
739 :
740 : /************************************************************************/
741 : /* SerializeToXML() */
742 : /************************************************************************/
743 :
744 4 : CPLXMLNode *VRTWarpedDataset::SerializeToXML( const char *pszVRTPath )
745 :
746 : {
747 : CPLXMLNode *psTree;
748 :
749 4 : psTree = VRTDataset::SerializeToXML( pszVRTPath );
750 :
751 4 : if( psTree == NULL )
752 0 : return psTree;
753 :
754 : /* -------------------------------------------------------------------- */
755 : /* Set subclass. */
756 : /* -------------------------------------------------------------------- */
757 : CPLCreateXMLNode(
758 : CPLCreateXMLNode( psTree, CXT_Attribute, "subClass" ),
759 4 : CXT_Text, "VRTWarpedDataset" );
760 :
761 : /* -------------------------------------------------------------------- */
762 : /* Serialize the block size. */
763 : /* -------------------------------------------------------------------- */
764 : CPLCreateXMLElementAndValue( psTree, "BlockXSize",
765 4 : CPLSPrintf( "%d", nBlockXSize ) );
766 : CPLCreateXMLElementAndValue( psTree, "BlockYSize",
767 4 : CPLSPrintf( "%d", nBlockYSize ) );
768 :
769 : /* -------------------------------------------------------------------- */
770 : /* Serialize the overview list. */
771 : /* -------------------------------------------------------------------- */
772 4 : if( nOverviewCount > 0 )
773 : {
774 : char *pszOverviewList;
775 : int iOverview;
776 :
777 1 : pszOverviewList = (char *) CPLMalloc(nOverviewCount*8 + 10);
778 1 : pszOverviewList[0] = '\0';
779 2 : for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
780 : {
781 : int nOvFactor;
782 :
783 : nOvFactor = (int)
784 : (0.5+GetRasterXSize()
785 1 : / (double) papoOverviews[iOverview]->GetRasterXSize());
786 :
787 : sprintf( pszOverviewList + strlen(pszOverviewList),
788 1 : "%d ", nOvFactor );
789 : }
790 :
791 1 : CPLCreateXMLElementAndValue( psTree, "OverviewList", pszOverviewList );
792 :
793 1 : CPLFree( pszOverviewList );
794 : }
795 :
796 : /* ==================================================================== */
797 : /* Serialize the warp options. */
798 : /* ==================================================================== */
799 : CPLXMLNode *psWOTree;
800 :
801 4 : if( poWarper != NULL )
802 : {
803 : /* -------------------------------------------------------------------- */
804 : /* We reset the destination dataset name so it doesn't get */
805 : /* written out in the serialize warp options. */
806 : /* -------------------------------------------------------------------- */
807 4 : char *pszSavedName = CPLStrdup(GetDescription());
808 4 : SetDescription("");
809 :
810 4 : psWOTree = GDALSerializeWarpOptions( poWarper->GetOptions() );
811 4 : CPLAddXMLChild( psTree, psWOTree );
812 :
813 4 : SetDescription( pszSavedName );
814 4 : CPLFree( pszSavedName );
815 :
816 : /* -------------------------------------------------------------------- */
817 : /* We need to consider making the source dataset relative to */
818 : /* the VRT file if possible. Adjust accordingly. */
819 : /* -------------------------------------------------------------------- */
820 4 : CPLXMLNode *psSDS = CPLGetXMLNode( psWOTree, "SourceDataset" );
821 : int bRelativeToVRT;
822 : char *pszRelativePath;
823 :
824 : pszRelativePath =
825 : CPLStrdup(
826 : CPLExtractRelativePath( pszVRTPath, psSDS->psChild->pszValue,
827 4 : &bRelativeToVRT ) );
828 :
829 4 : CPLFree( psSDS->psChild->pszValue );
830 4 : psSDS->psChild->pszValue = pszRelativePath;
831 :
832 : CPLCreateXMLNode(
833 : CPLCreateXMLNode( psSDS, CXT_Attribute, "relativeToVRT" ),
834 4 : CXT_Text, bRelativeToVRT ? "1" : "0" );
835 : }
836 :
837 4 : return psTree;
838 : }
839 :
840 : /************************************************************************/
841 : /* GetBlockSize() */
842 : /************************************************************************/
843 :
844 69 : void VRTWarpedDataset::GetBlockSize( int *pnBlockXSize, int *pnBlockYSize )
845 :
846 : {
847 69 : assert( NULL != pnBlockXSize );
848 69 : assert( NULL != pnBlockYSize );
849 :
850 69 : *pnBlockXSize = nBlockXSize;
851 69 : *pnBlockYSize = nBlockYSize;
852 69 : }
853 :
854 : /************************************************************************/
855 : /* ProcessBlock() */
856 : /* */
857 : /* Warp a single requested block, and then push each band of */
858 : /* the result into the block cache. */
859 : /************************************************************************/
860 :
861 65 : CPLErr VRTWarpedDataset::ProcessBlock( int iBlockX, int iBlockY )
862 :
863 : {
864 65 : if( poWarper == NULL )
865 0 : return CE_Failure;
866 :
867 65 : const GDALWarpOptions *psWO = poWarper->GetOptions();
868 :
869 : /* -------------------------------------------------------------------- */
870 : /* Allocate block of memory large enough to hold all the bands */
871 : /* for this block. */
872 : /* -------------------------------------------------------------------- */
873 : int iBand;
874 : GByte *pabyDstBuffer;
875 : int nDstBufferSize;
876 65 : int nWordSize = (GDALGetDataTypeSize(psWO->eWorkingDataType) / 8);
877 :
878 : // FIXME? : risk of overflow in multiplication if nBlockXSize or nBlockYSize are very large
879 65 : nDstBufferSize = nBlockXSize * nBlockYSize * psWO->nBandCount * nWordSize;
880 :
881 65 : pabyDstBuffer = (GByte *) VSIMalloc(nDstBufferSize);
882 :
883 65 : if( pabyDstBuffer == NULL )
884 : {
885 : CPLError( CE_Failure, CPLE_OutOfMemory,
886 : "Out of memory allocating %d byte buffer in VRTWarpedDataset::ProcessBlock()",
887 0 : nDstBufferSize );
888 0 : return CE_Failure;
889 : }
890 :
891 65 : memset( pabyDstBuffer, 0, nDstBufferSize );
892 :
893 : /* -------------------------------------------------------------------- */
894 : /* Process INIT_DEST option to initialize the buffer prior to */
895 : /* warping into it. */
896 : /* NOTE:The following code is 99% similar in gdalwarpoperation.cpp and */
897 : /* vrtwarped.cpp. Be careful to keep it in sync ! */
898 : /* -------------------------------------------------------------------- */
899 : const char *pszInitDest = CSLFetchNameValue( psWO->papszWarpOptions,
900 65 : "INIT_DEST" );
901 :
902 65 : if( pszInitDest != NULL && !EQUAL(pszInitDest, "") )
903 : {
904 : char **papszInitValues =
905 15 : CSLTokenizeStringComplex( pszInitDest, ",", FALSE, FALSE );
906 15 : int nInitCount = CSLCount(papszInitValues);
907 :
908 46 : for( iBand = 0; iBand < psWO->nBandCount; iBand++ )
909 : {
910 : double adfInitRealImag[2];
911 : GByte *pBandData;
912 31 : int nBandSize = nBlockXSize * nBlockYSize * nWordSize;
913 31 : const char *pszBandInit = papszInitValues[MIN(iBand,nInitCount-1)];
914 :
915 32 : if( EQUAL(pszBandInit,"NO_DATA")
916 : && psWO->padfDstNoDataReal != NULL )
917 : {
918 1 : adfInitRealImag[0] = psWO->padfDstNoDataReal[iBand];
919 1 : adfInitRealImag[1] = psWO->padfDstNoDataImag[iBand];
920 : }
921 : else
922 : {
923 : CPLStringToComplex( pszBandInit,
924 30 : adfInitRealImag + 0, adfInitRealImag + 1);
925 : }
926 :
927 31 : pBandData = ((GByte *) pabyDstBuffer) + iBand * nBandSize;
928 :
929 31 : if( psWO->eWorkingDataType == GDT_Byte )
930 : memset( pBandData,
931 37 : MAX(0,MIN(255,(int)adfInitRealImag[0])),
932 55 : nBandSize);
933 26 : else if( adfInitRealImag[0] == 0.0 && adfInitRealImag[1] == 0 )
934 : {
935 13 : memset( pBandData, 0, nBandSize );
936 : }
937 0 : else if( adfInitRealImag[1] == 0.0 )
938 : {
939 : GDALCopyWords( &adfInitRealImag, GDT_Float64, 0,
940 : pBandData,psWO->eWorkingDataType,nWordSize,
941 0 : nBlockXSize * nBlockYSize );
942 : }
943 : else
944 : {
945 : GDALCopyWords( &adfInitRealImag, GDT_CFloat64, 0,
946 : pBandData,psWO->eWorkingDataType,nWordSize,
947 0 : nBlockXSize * nBlockYSize );
948 : }
949 : }
950 :
951 15 : CSLDestroy( papszInitValues );
952 : }
953 :
954 : /* -------------------------------------------------------------------- */
955 : /* Warp into this buffer. */
956 : /* -------------------------------------------------------------------- */
957 : CPLErr eErr;
958 :
959 : eErr =
960 : poWarper->WarpRegionToBuffer(
961 : iBlockX * nBlockXSize, iBlockY * nBlockYSize,
962 : nBlockXSize, nBlockYSize,
963 65 : pabyDstBuffer, psWO->eWorkingDataType );
964 :
965 65 : if( eErr != CE_None )
966 : {
967 0 : VSIFree( pabyDstBuffer );
968 0 : return eErr;
969 : }
970 :
971 : /* -------------------------------------------------------------------- */
972 : /* Copy out into cache blocks for each band. */
973 : /* -------------------------------------------------------------------- */
974 152 : for( iBand = 0; iBand < psWO->nBandCount; iBand++ )
975 : {
976 : GDALRasterBand *poBand;
977 : GDALRasterBlock *poBlock;
978 :
979 87 : poBand = GetRasterBand(iBand+1);
980 87 : poBlock = poBand->GetLockedBlockRef( iBlockX, iBlockY, TRUE );
981 :
982 : CPLAssert( poBlock != NULL && poBlock->GetDataRef() != NULL );
983 :
984 : GDALCopyWords( pabyDstBuffer + iBand*nBlockXSize*nBlockYSize*nWordSize,
985 : psWO->eWorkingDataType, nWordSize,
986 : poBlock->GetDataRef(),
987 : poBlock->GetDataType(),
988 : GDALGetDataTypeSize(poBlock->GetDataType())/8,
989 87 : nBlockXSize * nBlockYSize );
990 :
991 87 : poBlock->DropLock();
992 : }
993 :
994 65 : VSIFree( pabyDstBuffer );
995 :
996 65 : return CE_None;
997 : }
998 :
999 : /************************************************************************/
1000 : /* AddBand() */
1001 : /************************************************************************/
1002 :
1003 5 : CPLErr VRTWarpedDataset::AddBand( GDALDataType eType, char **papszOptions )
1004 :
1005 : {
1006 : UNREFERENCED_PARAM( papszOptions );
1007 :
1008 : SetBand( GetRasterCount() + 1,
1009 5 : new VRTWarpedRasterBand( this, GetRasterCount() + 1, eType ) );
1010 :
1011 5 : return CE_None;
1012 : }
1013 :
1014 : /************************************************************************/
1015 : /* ==================================================================== */
1016 : /* VRTWarpedRasterBand */
1017 : /* ==================================================================== */
1018 : /************************************************************************/
1019 :
1020 : /************************************************************************/
1021 : /* VRTWarpedRasterBand() */
1022 : /************************************************************************/
1023 :
1024 69 : VRTWarpedRasterBand::VRTWarpedRasterBand( GDALDataset *poDS, int nBand,
1025 69 : GDALDataType eType )
1026 :
1027 : {
1028 69 : Initialize( poDS->GetRasterXSize(), poDS->GetRasterYSize() );
1029 :
1030 69 : this->poDS = poDS;
1031 69 : this->nBand = nBand;
1032 69 : this->eAccess = GA_Update;
1033 :
1034 : ((VRTWarpedDataset *) poDS)->GetBlockSize( &nBlockXSize,
1035 69 : &nBlockYSize );
1036 :
1037 69 : if( eType != GDT_Unknown )
1038 7 : this->eDataType = eType;
1039 69 : }
1040 :
1041 : /************************************************************************/
1042 : /* ~VRTWarpedRasterBand() */
1043 : /************************************************************************/
1044 :
1045 138 : VRTWarpedRasterBand::~VRTWarpedRasterBand()
1046 :
1047 : {
1048 69 : FlushCache();
1049 138 : }
1050 :
1051 : /************************************************************************/
1052 : /* IReadBlock() */
1053 : /************************************************************************/
1054 :
1055 65 : CPLErr VRTWarpedRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
1056 : void * pImage )
1057 :
1058 : {
1059 : CPLErr eErr;
1060 65 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1061 : GDALRasterBlock *poBlock;
1062 :
1063 65 : poBlock = GetLockedBlockRef( nBlockXOff, nBlockYOff, TRUE );
1064 :
1065 65 : eErr = poWDS->ProcessBlock( nBlockXOff, nBlockYOff );
1066 :
1067 65 : if( eErr == CE_None && pImage != poBlock->GetDataRef() )
1068 : {
1069 : int nDataBytes;
1070 : nDataBytes = (GDALGetDataTypeSize(poBlock->GetDataType()) / 8)
1071 0 : * poBlock->GetXSize() * poBlock->GetYSize();
1072 0 : memcpy( pImage, poBlock->GetDataRef(), nDataBytes );
1073 : }
1074 :
1075 65 : poBlock->DropLock();
1076 :
1077 65 : return eErr;
1078 : }
1079 :
1080 : /************************************************************************/
1081 : /* IWriteBlock() */
1082 : /************************************************************************/
1083 :
1084 8 : CPLErr VRTWarpedRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
1085 : void * pImage )
1086 :
1087 : {
1088 : CPLErr eErr;
1089 8 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1090 :
1091 : /* This is a bit tricky. In the case we are warping a VRTWarpedDataset */
1092 : /* with a destination alpha band, IWriteBlock can be called on that alpha */
1093 : /* band by GDALWarpDstAlphaMasker */
1094 : /* We don't need to do anything since the data will be kept in the block */
1095 : /* cache by VRTWarpedRasterBand::IReadBlock */
1096 8 : if (poWDS->poWarper->GetOptions()->nDstAlphaBand == nBand)
1097 : {
1098 8 : eErr = CE_None;
1099 : }
1100 : else
1101 : {
1102 : /* Otherwise, call the superclass method, that will fail of course */
1103 0 : eErr = VRTRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
1104 : }
1105 :
1106 8 : return eErr;
1107 : }
1108 :
1109 : /************************************************************************/
1110 : /* XMLInit() */
1111 : /************************************************************************/
1112 :
1113 62 : CPLErr VRTWarpedRasterBand::XMLInit( CPLXMLNode * psTree,
1114 : const char *pszVRTPath )
1115 :
1116 : {
1117 62 : return VRTRasterBand::XMLInit( psTree, pszVRTPath );
1118 : }
1119 :
1120 : /************************************************************************/
1121 : /* SerializeToXML() */
1122 : /************************************************************************/
1123 :
1124 6 : CPLXMLNode *VRTWarpedRasterBand::SerializeToXML( const char *pszVRTPath )
1125 :
1126 : {
1127 6 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML( pszVRTPath );
1128 :
1129 : /* -------------------------------------------------------------------- */
1130 : /* Set subclass. */
1131 : /* -------------------------------------------------------------------- */
1132 : CPLCreateXMLNode(
1133 : CPLCreateXMLNode( psTree, CXT_Attribute, "subClass" ),
1134 6 : CXT_Text, "VRTWarpedRasterBand" );
1135 :
1136 6 : return psTree;
1137 : }
1138 :
1139 : /************************************************************************/
1140 : /* GetOverviewCount() */
1141 : /************************************************************************/
1142 :
1143 1 : int VRTWarpedRasterBand::GetOverviewCount()
1144 :
1145 : {
1146 1 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1147 :
1148 1 : return poWDS->nOverviewCount;
1149 : }
1150 :
1151 : /************************************************************************/
1152 : /* GetOverview() */
1153 : /************************************************************************/
1154 :
1155 1 : GDALRasterBand *VRTWarpedRasterBand::GetOverview( int iOverview )
1156 :
1157 : {
1158 1 : VRTWarpedDataset *poWDS = (VRTWarpedDataset *) poDS;
1159 :
1160 1 : if( iOverview < 0 || iOverview >= poWDS->nOverviewCount )
1161 0 : return NULL;
1162 : else
1163 1 : return poWDS->papoOverviews[iOverview]->GetRasterBand( nBand );
1164 : }
1165 :
|
