Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL algorithms
4 : * Purpose: Apply vertical shift grid
5 : * Author: Even Rouault, even.rouault at spatialys.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2017, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 :
29 : #include "cpl_string.h"
30 : #include "gdal.h"
31 : #include "gdal_alg.h"
32 : #include "gdal_alg_priv.h"
33 : #include "gdal_priv.h"
34 : #include "gdal_utils.h"
35 : #include "gdalwarper.h"
36 : #include "vrtdataset.h"
37 : #include "ogr_spatialref.h"
38 :
39 : #include "proj.h"
40 :
41 : #include <limits>
42 :
43 : /************************************************************************/
44 : /* GDALApplyVSGDataset */
45 : /************************************************************************/
46 :
47 : class GDALApplyVSGDataset final : public GDALDataset
48 : {
49 : friend class GDALApplyVSGRasterBand;
50 :
51 : GDALDataset *m_poSrcDataset = nullptr;
52 : GDALDataset *m_poReprojectedGrid = nullptr;
53 : bool m_bInverse = false;
54 : double m_dfSrcUnitToMeter = 0.0;
55 : double m_dfDstUnitToMeter = 0.0;
56 :
57 : CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGDataset)
58 :
59 : public:
60 : GDALApplyVSGDataset(GDALDataset *poSrcDataset,
61 : GDALDataset *poReprojectedGrid, GDALDataType eDT,
62 : bool bInverse, double dfSrcUnitToMeter,
63 : double dfDstUnitToMeter, int nBlockSize);
64 : virtual ~GDALApplyVSGDataset();
65 :
66 : virtual int CloseDependentDatasets() override;
67 :
68 : virtual CPLErr GetGeoTransform(double *padfGeoTransform) override;
69 : virtual const OGRSpatialReference *GetSpatialRef() const override;
70 :
71 : bool IsInitOK();
72 : };
73 :
74 : /************************************************************************/
75 : /* GDALApplyVSGRasterBand */
76 : /************************************************************************/
77 :
78 : class GDALApplyVSGRasterBand final : public GDALRasterBand
79 : {
80 : friend class GDALApplyVSGDataset;
81 :
82 : float *m_pafSrcData = nullptr;
83 : float *m_pafGridData = nullptr;
84 :
85 : CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGRasterBand)
86 :
87 : public:
88 : GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize);
89 : virtual ~GDALApplyVSGRasterBand();
90 :
91 : virtual CPLErr IReadBlock(int nBlockXOff, int nBlockYOff,
92 : void *pData) override;
93 : virtual double GetNoDataValue(int *pbSuccess) override;
94 : };
95 :
96 : /************************************************************************/
97 : /* GDALApplyVSGDataset() */
98 : /************************************************************************/
99 :
100 13 : GDALApplyVSGDataset::GDALApplyVSGDataset(GDALDataset *poSrcDataset,
101 : GDALDataset *poReprojectedGrid,
102 : GDALDataType eDT, bool bInverse,
103 : double dfSrcUnitToMeter,
104 : double dfDstUnitToMeter,
105 13 : int nBlockSize)
106 : : m_poSrcDataset(poSrcDataset), m_poReprojectedGrid(poReprojectedGrid),
107 : m_bInverse(bInverse), m_dfSrcUnitToMeter(dfSrcUnitToMeter),
108 13 : m_dfDstUnitToMeter(dfDstUnitToMeter)
109 : {
110 13 : m_poSrcDataset->Reference();
111 13 : m_poReprojectedGrid->Reference();
112 :
113 13 : nRasterXSize = poSrcDataset->GetRasterXSize();
114 13 : nRasterYSize = poSrcDataset->GetRasterYSize();
115 13 : SetBand(1, new GDALApplyVSGRasterBand(eDT, nBlockSize));
116 13 : }
117 :
118 : /************************************************************************/
119 : /* ~GDALApplyVSGDataset() */
120 : /************************************************************************/
121 :
122 26 : GDALApplyVSGDataset::~GDALApplyVSGDataset()
123 : {
124 13 : GDALApplyVSGDataset::CloseDependentDatasets();
125 26 : }
126 :
127 : /************************************************************************/
128 : /* CloseDependentDatasets() */
129 : /************************************************************************/
130 :
131 13 : int GDALApplyVSGDataset::CloseDependentDatasets()
132 : {
133 13 : bool bRet = false;
134 13 : if (m_poSrcDataset != nullptr)
135 : {
136 13 : if (m_poSrcDataset->ReleaseRef())
137 : {
138 8 : bRet = true;
139 : }
140 13 : m_poSrcDataset = nullptr;
141 : }
142 13 : if (m_poReprojectedGrid != nullptr)
143 : {
144 13 : if (m_poReprojectedGrid->ReleaseRef())
145 : {
146 13 : bRet = true;
147 : }
148 13 : m_poReprojectedGrid = nullptr;
149 : }
150 13 : return bRet;
151 : }
152 :
153 : /************************************************************************/
154 : /* GetGeoTransform() */
155 : /************************************************************************/
156 :
157 2 : CPLErr GDALApplyVSGDataset::GetGeoTransform(double *padfGeoTransform)
158 : {
159 2 : return m_poSrcDataset->GetGeoTransform(padfGeoTransform);
160 : }
161 :
162 : /************************************************************************/
163 : /* GetSpatialRef() */
164 : /************************************************************************/
165 :
166 2 : const OGRSpatialReference *GDALApplyVSGDataset::GetSpatialRef() const
167 : {
168 2 : return m_poSrcDataset->GetSpatialRef();
169 : }
170 :
171 : /************************************************************************/
172 : /* IsInitOK() */
173 : /************************************************************************/
174 :
175 13 : bool GDALApplyVSGDataset::IsInitOK()
176 : {
177 : GDALApplyVSGRasterBand *poBand =
178 13 : reinterpret_cast<GDALApplyVSGRasterBand *>(GetRasterBand(1));
179 13 : return poBand->m_pafSrcData != nullptr && poBand->m_pafGridData != nullptr;
180 : }
181 :
182 : /************************************************************************/
183 : /* GDALApplyVSGRasterBand() */
184 : /************************************************************************/
185 :
186 13 : GDALApplyVSGRasterBand::GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize)
187 : {
188 13 : eDataType = eDT;
189 13 : nBlockXSize = nBlockSize;
190 13 : nBlockYSize = nBlockSize;
191 13 : m_pafSrcData = static_cast<float *>(
192 13 : VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));
193 13 : m_pafGridData = static_cast<float *>(
194 13 : VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));
195 13 : }
196 :
197 : /************************************************************************/
198 : /* ~GDALApplyVSGRasterBand() */
199 : /************************************************************************/
200 :
201 26 : GDALApplyVSGRasterBand::~GDALApplyVSGRasterBand()
202 : {
203 13 : VSIFree(m_pafSrcData);
204 13 : VSIFree(m_pafGridData);
205 26 : }
206 :
207 : /************************************************************************/
208 : /* GetNoDataValue() */
209 : /************************************************************************/
210 :
211 19 : double GDALApplyVSGRasterBand::GetNoDataValue(int *pbSuccess)
212 : {
213 19 : GDALApplyVSGDataset *poGDS = reinterpret_cast<GDALApplyVSGDataset *>(poDS);
214 19 : return poGDS->m_poSrcDataset->GetRasterBand(1)->GetNoDataValue(pbSuccess);
215 : }
216 :
217 : /************************************************************************/
218 : /* IReadBlock() */
219 : /************************************************************************/
220 :
221 17 : CPLErr GDALApplyVSGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
222 : void *pData)
223 : {
224 17 : GDALApplyVSGDataset *poGDS = reinterpret_cast<GDALApplyVSGDataset *>(poDS);
225 :
226 17 : const int nXOff = nBlockXOff * nBlockXSize;
227 34 : const int nReqXSize = (nXOff > nRasterXSize - nBlockXSize)
228 17 : ? nRasterXSize - nXOff
229 : : nBlockXSize;
230 17 : const int nYOff = nBlockYOff * nBlockYSize;
231 34 : const int nReqYSize = (nYOff > nRasterYSize - nBlockYSize)
232 17 : ? nRasterYSize - nYOff
233 : : nBlockYSize;
234 :
235 17 : CPLErr eErr = poGDS->m_poSrcDataset->GetRasterBand(1)->RasterIO(
236 17 : GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafSrcData, nReqXSize,
237 17 : nReqYSize, GDT_Float32, sizeof(float), nBlockXSize * sizeof(float),
238 : nullptr);
239 17 : if (eErr == CE_None)
240 34 : eErr = poGDS->m_poReprojectedGrid->GetRasterBand(1)->RasterIO(
241 17 : GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafGridData,
242 : nReqXSize, nReqYSize, GDT_Float32, sizeof(float),
243 17 : nBlockXSize * sizeof(float), nullptr);
244 17 : if (eErr == CE_None)
245 : {
246 17 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
247 17 : int bHasNoData = FALSE;
248 17 : float fNoDataValue = static_cast<float>(GetNoDataValue(&bHasNoData));
249 279 : for (int iY = 0; iY < nReqYSize; iY++)
250 : {
251 4666 : for (int iX = 0; iX < nReqXSize; iX++)
252 : {
253 4404 : const float fSrcVal = m_pafSrcData[iY * nBlockXSize + iX];
254 4404 : const float fGridVal = m_pafGridData[iY * nBlockXSize + iX];
255 4404 : if (bHasNoData && fSrcVal == fNoDataValue)
256 : {
257 : }
258 4403 : else if (CPLIsInf(fGridVal))
259 : {
260 2 : CPLError(CE_Failure, CPLE_AppDefined,
261 : "Missing vertical grid value at source (%d,%d)",
262 : nXOff + iX, nYOff + iY);
263 2 : return CE_Failure;
264 : }
265 4401 : else if (poGDS->m_bInverse)
266 : {
267 800 : m_pafSrcData[iY * nBlockXSize + iX] = static_cast<float>(
268 800 : (fSrcVal * poGDS->m_dfSrcUnitToMeter - fGridVal) /
269 800 : poGDS->m_dfDstUnitToMeter);
270 : }
271 : else
272 : {
273 3601 : m_pafSrcData[iY * nBlockXSize + iX] = static_cast<float>(
274 3601 : (fSrcVal * poGDS->m_dfSrcUnitToMeter + fGridVal) /
275 3601 : poGDS->m_dfDstUnitToMeter);
276 : }
277 : }
278 262 : GDALCopyWords(
279 262 : m_pafSrcData + iY * nBlockXSize, GDT_Float32, sizeof(float),
280 262 : static_cast<GByte *>(pData) + iY * nBlockXSize * nDTSize,
281 : eDataType, nDTSize, nReqXSize);
282 : }
283 : }
284 :
285 15 : return eErr;
286 : }
287 :
288 : /************************************************************************/
289 : /* GDALApplyVerticalShiftGrid() */
290 : /************************************************************************/
291 :
292 : /** Apply a vertical shift grid to a source (DEM typically) dataset.
293 : *
294 : * hGridDataset will typically use WGS84 as horizontal datum (but this is
295 : * not a requirement) and its values are the values to add to go from geoid
296 : * elevations to WGS84 ellipsoidal heights.
297 : *
298 : * hGridDataset will be on-the-fly reprojected and resampled to the projection
299 : * and resolution of hSrcDataset, using bilinear resampling by default.
300 : *
301 : * Both hSrcDataset and hGridDataset must be single band datasets, and have
302 : * a valid geotransform and projection.
303 : *
304 : * On success, a reference will be taken on hSrcDataset and hGridDataset.
305 : * Reference counting semantics on the source and grid datasets should be
306 : * honoured. That is, don't just GDALClose() it, unless it was opened with
307 : * GDALOpenShared(), but rather use GDALReleaseDataset() if wanting to
308 : * immediately release the reference(s) and make the returned dataset the
309 : * owner of them.
310 : *
311 : * Valid use cases:
312 : *
313 : * \code
314 : * hSrcDataset = GDALOpen(...)
315 : * hGridDataset = GDALOpen(...)
316 : * hDstDataset = GDALApplyVerticalShiftGrid(hSrcDataset, hGridDataset, ...)
317 : * GDALReleaseDataset(hSrcDataset);
318 : * GDALReleaseDataset(hGridDataset);
319 : * if( hDstDataset )
320 : * {
321 : * // Do things with hDstDataset
322 : * GDALClose(hDstDataset) // will close hSrcDataset and hGridDataset
323 : * }
324 : * \endcode
325 :
326 : *
327 : * @param hSrcDataset source (DEM) dataset. Must not be NULL.
328 : * @param hGridDataset vertical grid shift dataset. Must not be NULL.
329 : * @param bInverse if set to FALSE, hGridDataset values will be added to
330 : * hSrcDataset. If set to TRUE, they will be subtracted.
331 : * @param dfSrcUnitToMeter the factor to convert values from hSrcDataset to
332 : * meters (1.0 if source values are in meter).
333 : * @param dfDstUnitToMeter the factor to convert shifted values from meter
334 : * (1.0 if output values must be in meter).
335 : * @param papszOptions list of options, or NULL. Supported options are:
336 : * <ul>
337 : * <li>RESAMPLING=NEAREST/BILINEAR/CUBIC. Defaults to BILINEAR.</li>
338 : * <li>MAX_ERROR=val. Maximum error measured in input pixels that is allowed in
339 : * approximating the transformation (0.0 for exact calculations). Defaults
340 : * to 0.125</li>
341 : * <li>DATATYPE=Byte/UInt16/Int16/Float32/Float64. Output data type. If not
342 : * specified will be the same as the one of hSrcDataset.
343 : * <li>ERROR_ON_MISSING_VERT_SHIFT=YES/NO. Whether a missing/nodata value in
344 : * hGridDataset should cause I/O requests to fail. Default is NO (in which case
345 : * 0 will be used)
346 : * <li>SRC_SRS=srs_def. Override projection on hSrcDataset;
347 : * </ul>
348 : *
349 : * @return a new dataset corresponding to hSrcDataset adjusted with
350 : * hGridDataset, or NULL. If not NULL, it must be closed with GDALClose().
351 : *
352 : * @since GDAL 2.2
353 : * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used
354 : * by gdalwarp initially, but is no longer needed.
355 : */
356 23 : GDALDatasetH GDALApplyVerticalShiftGrid(GDALDatasetH hSrcDataset,
357 : GDALDatasetH hGridDataset, int bInverse,
358 : double dfSrcUnitToMeter,
359 : double dfDstUnitToMeter,
360 : const char *const *papszOptions)
361 : {
362 23 : VALIDATE_POINTER1(hSrcDataset, "GDALApplyVerticalShiftGrid", nullptr);
363 23 : VALIDATE_POINTER1(hGridDataset, "GDALApplyVerticalShiftGrid", nullptr);
364 :
365 : double adfSrcGT[6];
366 23 : if (GDALGetGeoTransform(hSrcDataset, adfSrcGT) != CE_None)
367 : {
368 1 : CPLError(CE_Failure, CPLE_NotSupported,
369 : "Source dataset has no geotransform.");
370 1 : return nullptr;
371 : }
372 22 : const char *pszSrcProjection = CSLFetchNameValue(papszOptions, "SRC_SRS");
373 44 : OGRSpatialReference oSrcSRS;
374 22 : if (pszSrcProjection != nullptr && pszSrcProjection[0] != '\0')
375 : {
376 0 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
377 0 : oSrcSRS.SetFromUserInput(pszSrcProjection);
378 : }
379 : else
380 : {
381 22 : auto poSRS = GDALDataset::FromHandle(hSrcDataset)->GetSpatialRef();
382 22 : if (poSRS)
383 21 : oSrcSRS = *poSRS;
384 : }
385 :
386 22 : if (oSrcSRS.IsCompound())
387 : {
388 0 : oSrcSRS.StripVertical();
389 : }
390 :
391 22 : if (oSrcSRS.IsEmpty())
392 : {
393 1 : CPLError(CE_Failure, CPLE_NotSupported,
394 : "Source dataset has no projection.");
395 1 : return nullptr;
396 : }
397 21 : if (GDALGetRasterCount(hSrcDataset) != 1)
398 : {
399 1 : CPLError(CE_Failure, CPLE_NotSupported,
400 : "Only single band source dataset is supported.");
401 1 : return nullptr;
402 : }
403 :
404 : double adfGridGT[6];
405 20 : if (GDALGetGeoTransform(hGridDataset, adfGridGT) != CE_None)
406 : {
407 1 : CPLError(CE_Failure, CPLE_NotSupported,
408 : "Grid dataset has no geotransform.");
409 1 : return nullptr;
410 : }
411 :
412 19 : OGRSpatialReferenceH hGridSRS = GDALGetSpatialRef(hGridDataset);
413 19 : if (hGridSRS == nullptr)
414 : {
415 1 : CPLError(CE_Failure, CPLE_NotSupported,
416 : "Grid dataset has no projection.");
417 1 : return nullptr;
418 : }
419 18 : if (GDALGetRasterCount(hGridDataset) != 1)
420 : {
421 1 : CPLError(CE_Failure, CPLE_NotSupported,
422 : "Only single band grid dataset is supported.");
423 1 : return nullptr;
424 : }
425 :
426 17 : GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDataset, 1));
427 17 : const char *pszDataType = CSLFetchNameValue(papszOptions, "DATATYPE");
428 17 : if (pszDataType)
429 2 : eDT = GDALGetDataTypeByName(pszDataType);
430 17 : if (eDT == GDT_Unknown)
431 : {
432 1 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid DATATYPE=%s",
433 : pszDataType);
434 1 : return nullptr;
435 : }
436 :
437 16 : const int nSrcXSize = GDALGetRasterXSize(hSrcDataset);
438 16 : const int nSrcYSize = GDALGetRasterYSize(hSrcDataset);
439 :
440 16 : double dfWestLongitudeDeg = 0.0;
441 16 : double dfSouthLatitudeDeg = 0.0;
442 16 : double dfEastLongitudeDeg = 0.0;
443 16 : double dfNorthLatitudeDeg = 0.0;
444 16 : GDALComputeAreaOfInterest(&oSrcSRS, adfSrcGT, nSrcXSize, nSrcYSize,
445 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
446 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
447 :
448 32 : CPLStringList aosOptions;
449 16 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
450 1 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
451 : {
452 : aosOptions.SetNameValue(
453 : "AREA_OF_INTEREST",
454 : CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
455 : dfSouthLatitudeDeg, dfEastLongitudeDeg,
456 15 : dfNorthLatitudeDeg));
457 : }
458 16 : void *hTransform = GDALCreateGenImgProjTransformer4(
459 : hGridSRS, adfGridGT, OGRSpatialReference::ToHandle(&oSrcSRS), adfSrcGT,
460 16 : aosOptions.List());
461 16 : if (hTransform == nullptr)
462 3 : return nullptr;
463 13 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
464 13 : psWO->hSrcDS = hGridDataset;
465 13 : psWO->eResampleAlg = GRA_Bilinear;
466 13 : const char *pszResampling = CSLFetchNameValue(papszOptions, "RESAMPLING");
467 13 : if (pszResampling)
468 : {
469 3 : if (EQUAL(pszResampling, "NEAREST"))
470 1 : psWO->eResampleAlg = GRA_NearestNeighbour;
471 2 : else if (EQUAL(pszResampling, "BILINEAR"))
472 1 : psWO->eResampleAlg = GRA_Bilinear;
473 1 : else if (EQUAL(pszResampling, "CUBIC"))
474 1 : psWO->eResampleAlg = GRA_Cubic;
475 : }
476 13 : psWO->eWorkingDataType = GDT_Float32;
477 13 : int bHasNoData = FALSE;
478 13 : const double dfSrcNoData = GDALGetRasterNoDataValue(
479 : GDALGetRasterBand(hGridDataset, 1), &bHasNoData);
480 13 : if (bHasNoData)
481 : {
482 2 : psWO->padfSrcNoDataReal =
483 2 : static_cast<double *>(CPLMalloc(sizeof(double)));
484 2 : psWO->padfSrcNoDataReal[0] = dfSrcNoData;
485 : }
486 :
487 13 : psWO->padfDstNoDataReal = static_cast<double *>(CPLMalloc(sizeof(double)));
488 : const bool bErrorOnMissingShift =
489 13 : CPLFetchBool(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false);
490 13 : psWO->padfDstNoDataReal[0] =
491 13 : (bErrorOnMissingShift) ? -std::numeric_limits<float>::infinity() : 0.0;
492 13 : psWO->papszWarpOptions =
493 13 : CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
494 :
495 13 : psWO->pfnTransformer = GDALGenImgProjTransform;
496 13 : psWO->pTransformerArg = hTransform;
497 : const double dfMaxError =
498 13 : CPLAtof(CSLFetchNameValueDef(papszOptions, "MAX_ERROR", "0.125"));
499 13 : if (dfMaxError > 0.0)
500 : {
501 13 : psWO->pTransformerArg = GDALCreateApproxTransformer(
502 : psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);
503 13 : psWO->pfnTransformer = GDALApproxTransform;
504 13 : GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);
505 : }
506 13 : psWO->nBandCount = 1;
507 13 : psWO->panSrcBands = static_cast<int *>(CPLMalloc(sizeof(int)));
508 13 : psWO->panSrcBands[0] = 1;
509 13 : psWO->panDstBands = static_cast<int *>(CPLMalloc(sizeof(int)));
510 13 : psWO->panDstBands[0] = 1;
511 :
512 : VRTWarpedDataset *poReprojectedGrid =
513 13 : new VRTWarpedDataset(nSrcXSize, nSrcYSize);
514 : // This takes a reference on hGridDataset
515 13 : CPLErr eErr = poReprojectedGrid->Initialize(psWO);
516 13 : CPLAssert(eErr == CE_None);
517 13 : CPL_IGNORE_RET_VAL(eErr);
518 13 : GDALDestroyWarpOptions(psWO);
519 13 : poReprojectedGrid->SetGeoTransform(adfSrcGT);
520 13 : poReprojectedGrid->AddBand(GDT_Float32, nullptr);
521 :
522 : GDALApplyVSGDataset *poOutDS = new GDALApplyVSGDataset(
523 13 : GDALDataset::FromHandle(hSrcDataset), poReprojectedGrid, eDT,
524 13 : CPL_TO_BOOL(bInverse), dfSrcUnitToMeter, dfDstUnitToMeter,
525 : // Undocumented option. For testing only
526 13 : atoi(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", "256")));
527 :
528 13 : poReprojectedGrid->ReleaseRef();
529 :
530 13 : if (!poOutDS->IsInitOK())
531 : {
532 1 : delete poOutDS;
533 1 : return nullptr;
534 : }
535 12 : poOutDS->SetDescription(GDALGetDescription(hSrcDataset));
536 12 : return reinterpret_cast<GDALDatasetH>(poOutDS);
537 : }
538 :
539 : /************************************************************************/
540 : /* GetProj4Filename() */
541 : /************************************************************************/
542 :
543 0 : static CPLString GetProj4Filename(const char *pszFilename)
544 : {
545 0 : CPLString osFilename;
546 :
547 : /* or fixed path: /name, ./name or ../name */
548 0 : if (!CPLIsFilenameRelative(pszFilename) || *pszFilename == '.')
549 : {
550 0 : return pszFilename;
551 : }
552 :
553 0 : PJ_GRID_INFO info = proj_grid_info(pszFilename);
554 0 : if (info.filename[0])
555 : {
556 0 : osFilename = info.filename;
557 : }
558 :
559 0 : return osFilename;
560 : }
561 :
562 : /************************************************************************/
563 : /* GDALOpenVerticalShiftGrid() */
564 : /************************************************************************/
565 :
566 : /** Load proj.4 geoidgrids as GDAL dataset
567 : *
568 : * @param pszProj4Geoidgrids Value of proj.4 geoidgrids parameter.
569 : * @param pbError If not NULL, the pointed value will be set to TRUE if an
570 : * error occurred.
571 : *
572 : * @return a dataset. If not NULL, it must be closed with GDALClose().
573 : *
574 : * @since GDAL 2.2
575 : * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used
576 : * by gdalwarp initially, but is no longer needed.
577 : */
578 0 : GDALDatasetH GDALOpenVerticalShiftGrid(const char *pszProj4Geoidgrids,
579 : int *pbError)
580 : {
581 0 : char **papszGrids = CSLTokenizeString2(pszProj4Geoidgrids, ",", 0);
582 0 : const int nGridCount = CSLCount(papszGrids);
583 0 : if (nGridCount == 1)
584 : {
585 0 : CSLDestroy(papszGrids);
586 :
587 0 : bool bMissingOk = false;
588 0 : if (*pszProj4Geoidgrids == '@')
589 : {
590 0 : pszProj4Geoidgrids++;
591 0 : bMissingOk = true;
592 : }
593 0 : const CPLString osFilename(GetProj4Filename(pszProj4Geoidgrids));
594 0 : const char *const papszOpenOptions[] = {
595 : "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES", nullptr};
596 : GDALDatasetH hDS =
597 0 : GDALOpenEx(osFilename, 0, nullptr, papszOpenOptions, nullptr);
598 0 : if (hDS == nullptr)
599 : {
600 0 : CPLDebug("GDAL", "Cannot find file corresponding to %s",
601 : pszProj4Geoidgrids);
602 : }
603 0 : if (pbError)
604 0 : *pbError = (!bMissingOk && hDS == nullptr);
605 0 : return hDS;
606 : }
607 :
608 0 : CPLStringList aosFilenames;
609 0 : for (int i = nGridCount - 1; i >= 0; i--)
610 : {
611 0 : const char *pszName = papszGrids[i];
612 0 : bool bMissingOk = false;
613 0 : if (*pszName == '@')
614 : {
615 0 : pszName++;
616 0 : bMissingOk = true;
617 : }
618 0 : const CPLString osFilename(GetProj4Filename(pszName));
619 : VSIStatBufL sStat;
620 0 : if (osFilename.empty() || VSIStatL(osFilename, &sStat) != 0)
621 : {
622 0 : CPLDebug("GDAL", "Cannot find file corresponding to %s", pszName);
623 0 : if (!bMissingOk)
624 : {
625 0 : if (pbError)
626 0 : *pbError = true;
627 0 : CSLDestroy(papszGrids);
628 0 : return nullptr;
629 : }
630 : }
631 : else
632 : {
633 0 : aosFilenames.AddString(osFilename);
634 : }
635 : }
636 :
637 0 : CSLDestroy(papszGrids);
638 :
639 0 : if (aosFilenames.empty())
640 : {
641 0 : if (pbError)
642 0 : *pbError = false;
643 0 : return nullptr;
644 : }
645 :
646 0 : char **papszArgv = nullptr;
647 0 : papszArgv = CSLAddString(papszArgv, "-resolution");
648 0 : papszArgv = CSLAddString(papszArgv, "highest");
649 0 : papszArgv = CSLAddString(papszArgv, "-vrtnodata");
650 0 : papszArgv = CSLAddString(papszArgv, "-inf");
651 0 : papszArgv = CSLAddString(papszArgv, "-oo");
652 : papszArgv =
653 0 : CSLAddString(papszArgv, "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES");
654 0 : GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(papszArgv, nullptr);
655 0 : CSLDestroy(papszArgv);
656 0 : GDALDatasetH hDS = GDALBuildVRT("", aosFilenames.size(), nullptr,
657 0 : aosFilenames.List(), psOptions, nullptr);
658 0 : GDALBuildVRTOptionsFree(psOptions);
659 0 : if (pbError)
660 0 : *pbError = hDS != nullptr;
661 0 : return hDS;
662 : }
|
