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