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 : * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used
337 : * by gdalwarp initially, but is no longer needed.
338 : */
339 23 : GDALDatasetH GDALApplyVerticalShiftGrid(GDALDatasetH hSrcDataset,
340 : GDALDatasetH hGridDataset, int bInverse,
341 : double dfSrcUnitToMeter,
342 : double dfDstUnitToMeter,
343 : const char *const *papszOptions)
344 : {
345 23 : VALIDATE_POINTER1(hSrcDataset, "GDALApplyVerticalShiftGrid", nullptr);
346 23 : VALIDATE_POINTER1(hGridDataset, "GDALApplyVerticalShiftGrid", nullptr);
347 :
348 23 : GDALGeoTransform srcGT;
349 23 : if (GDALDataset::FromHandle(hSrcDataset)->GetGeoTransform(srcGT) != CE_None)
350 : {
351 1 : CPLError(CE_Failure, CPLE_NotSupported,
352 : "Source dataset has no geotransform.");
353 1 : return nullptr;
354 : }
355 22 : const char *pszSrcProjection = CSLFetchNameValue(papszOptions, "SRC_SRS");
356 44 : OGRSpatialReference oSrcSRS;
357 22 : if (pszSrcProjection != nullptr && pszSrcProjection[0] != '\0')
358 : {
359 0 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
360 0 : oSrcSRS.SetFromUserInput(pszSrcProjection);
361 : }
362 : else
363 : {
364 22 : auto poSRS = GDALDataset::FromHandle(hSrcDataset)->GetSpatialRef();
365 22 : if (poSRS)
366 21 : oSrcSRS = *poSRS;
367 : }
368 :
369 22 : if (oSrcSRS.IsCompound())
370 : {
371 0 : oSrcSRS.StripVertical();
372 : }
373 :
374 22 : if (oSrcSRS.IsEmpty())
375 : {
376 1 : CPLError(CE_Failure, CPLE_NotSupported,
377 : "Source dataset has no projection.");
378 1 : return nullptr;
379 : }
380 21 : if (GDALGetRasterCount(hSrcDataset) != 1)
381 : {
382 1 : CPLError(CE_Failure, CPLE_NotSupported,
383 : "Only single band source dataset is supported.");
384 1 : return nullptr;
385 : }
386 :
387 20 : GDALGeoTransform gridGT;
388 20 : if (GDALDataset::FromHandle(hGridDataset)->GetGeoTransform(gridGT) !=
389 : CE_None)
390 : {
391 1 : CPLError(CE_Failure, CPLE_NotSupported,
392 : "Grid dataset has no geotransform.");
393 1 : return nullptr;
394 : }
395 :
396 19 : OGRSpatialReferenceH hGridSRS = GDALGetSpatialRef(hGridDataset);
397 19 : if (hGridSRS == nullptr)
398 : {
399 1 : CPLError(CE_Failure, CPLE_NotSupported,
400 : "Grid dataset has no projection.");
401 1 : return nullptr;
402 : }
403 18 : if (GDALGetRasterCount(hGridDataset) != 1)
404 : {
405 1 : CPLError(CE_Failure, CPLE_NotSupported,
406 : "Only single band grid dataset is supported.");
407 1 : return nullptr;
408 : }
409 :
410 17 : GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDataset, 1));
411 17 : const char *pszDataType = CSLFetchNameValue(papszOptions, "DATATYPE");
412 17 : if (pszDataType)
413 2 : eDT = GDALGetDataTypeByName(pszDataType);
414 17 : if (eDT == GDT_Unknown)
415 : {
416 1 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid DATATYPE=%s",
417 : pszDataType);
418 1 : return nullptr;
419 : }
420 :
421 16 : const int nSrcXSize = GDALGetRasterXSize(hSrcDataset);
422 16 : const int nSrcYSize = GDALGetRasterYSize(hSrcDataset);
423 :
424 16 : double dfWestLongitudeDeg = 0.0;
425 16 : double dfSouthLatitudeDeg = 0.0;
426 16 : double dfEastLongitudeDeg = 0.0;
427 16 : double dfNorthLatitudeDeg = 0.0;
428 16 : GDALComputeAreaOfInterest(&oSrcSRS, srcGT.data(), nSrcXSize, nSrcYSize,
429 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
430 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
431 :
432 32 : CPLStringList aosOptions;
433 16 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
434 1 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
435 : {
436 : aosOptions.SetNameValue(
437 : "AREA_OF_INTEREST",
438 : CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
439 : dfSouthLatitudeDeg, dfEastLongitudeDeg,
440 15 : dfNorthLatitudeDeg));
441 : }
442 32 : void *hTransform = GDALCreateGenImgProjTransformer4(
443 16 : hGridSRS, gridGT.data(), OGRSpatialReference::ToHandle(&oSrcSRS),
444 16 : srcGT.data(), aosOptions.List());
445 16 : if (hTransform == nullptr)
446 3 : return nullptr;
447 13 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
448 13 : psWO->hSrcDS = hGridDataset;
449 13 : psWO->eResampleAlg = GRA_Bilinear;
450 13 : const char *pszResampling = CSLFetchNameValue(papszOptions, "RESAMPLING");
451 13 : if (pszResampling)
452 : {
453 3 : if (EQUAL(pszResampling, "NEAREST"))
454 1 : psWO->eResampleAlg = GRA_NearestNeighbour;
455 2 : else if (EQUAL(pszResampling, "BILINEAR"))
456 1 : psWO->eResampleAlg = GRA_Bilinear;
457 1 : else if (EQUAL(pszResampling, "CUBIC"))
458 1 : psWO->eResampleAlg = GRA_Cubic;
459 : }
460 13 : psWO->eWorkingDataType = GDT_Float32;
461 13 : int bHasNoData = FALSE;
462 13 : const double dfSrcNoData = GDALGetRasterNoDataValue(
463 : GDALGetRasterBand(hGridDataset, 1), &bHasNoData);
464 13 : if (bHasNoData)
465 : {
466 2 : psWO->padfSrcNoDataReal =
467 2 : static_cast<double *>(CPLMalloc(sizeof(double)));
468 2 : psWO->padfSrcNoDataReal[0] = dfSrcNoData;
469 : }
470 :
471 13 : psWO->padfDstNoDataReal = static_cast<double *>(CPLMalloc(sizeof(double)));
472 : const bool bErrorOnMissingShift =
473 13 : CPLFetchBool(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false);
474 13 : psWO->padfDstNoDataReal[0] =
475 : (bErrorOnMissingShift)
476 13 : ? static_cast<double>(-std::numeric_limits<float>::infinity())
477 : : 0.0;
478 13 : psWO->papszWarpOptions =
479 13 : CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
480 :
481 13 : psWO->pfnTransformer = GDALGenImgProjTransform;
482 13 : psWO->pTransformerArg = hTransform;
483 : const double dfMaxError =
484 13 : CPLAtof(CSLFetchNameValueDef(papszOptions, "MAX_ERROR", "0.125"));
485 13 : if (dfMaxError > 0.0)
486 : {
487 13 : psWO->pTransformerArg = GDALCreateApproxTransformer(
488 : psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);
489 13 : psWO->pfnTransformer = GDALApproxTransform;
490 13 : GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);
491 : }
492 13 : psWO->nBandCount = 1;
493 13 : psWO->panSrcBands = static_cast<int *>(CPLMalloc(sizeof(int)));
494 13 : psWO->panSrcBands[0] = 1;
495 13 : psWO->panDstBands = static_cast<int *>(CPLMalloc(sizeof(int)));
496 13 : psWO->panDstBands[0] = 1;
497 :
498 : VRTWarpedDataset *poReprojectedGrid =
499 13 : new VRTWarpedDataset(nSrcXSize, nSrcYSize);
500 : // This takes a reference on hGridDataset
501 13 : CPLErr eErr = poReprojectedGrid->Initialize(psWO);
502 13 : CPLAssert(eErr == CE_None);
503 13 : CPL_IGNORE_RET_VAL(eErr);
504 13 : GDALDestroyWarpOptions(psWO);
505 13 : poReprojectedGrid->SetGeoTransform(srcGT);
506 13 : poReprojectedGrid->AddBand(GDT_Float32, nullptr);
507 :
508 : GDALApplyVSGDataset *poOutDS = new GDALApplyVSGDataset(
509 13 : GDALDataset::FromHandle(hSrcDataset), poReprojectedGrid, eDT,
510 13 : CPL_TO_BOOL(bInverse), static_cast<float>(dfSrcUnitToMeter),
511 : static_cast<float>(dfDstUnitToMeter),
512 : // Undocumented option. For testing only
513 13 : atoi(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", "256")));
514 :
515 13 : poReprojectedGrid->ReleaseRef();
516 :
517 13 : if (!poOutDS->IsInitOK())
518 : {
519 1 : delete poOutDS;
520 1 : return nullptr;
521 : }
522 12 : poOutDS->SetDescription(GDALGetDescription(hSrcDataset));
523 12 : return reinterpret_cast<GDALDatasetH>(poOutDS);
524 : }
525 :
526 : /************************************************************************/
527 : /* GetProj4Filename() */
528 : /************************************************************************/
529 :
530 0 : static CPLString GetProj4Filename(const char *pszFilename)
531 : {
532 0 : CPLString osFilename;
533 :
534 : /* or fixed path: /name, ./name or ../name */
535 0 : if (!CPLIsFilenameRelative(pszFilename) || *pszFilename == '.')
536 : {
537 0 : return pszFilename;
538 : }
539 :
540 0 : PJ_GRID_INFO info = proj_grid_info(pszFilename);
541 0 : if (info.filename[0])
542 : {
543 0 : osFilename = info.filename;
544 : }
545 :
546 0 : return osFilename;
547 : }
548 :
549 : /************************************************************************/
550 : /* GDALOpenVerticalShiftGrid() */
551 : /************************************************************************/
552 :
553 : /** Load proj.4 geoidgrids as GDAL dataset
554 : *
555 : * @param pszProj4Geoidgrids Value of proj.4 geoidgrids parameter.
556 : * @param pbError If not NULL, the pointed value will be set to TRUE if an
557 : * error occurred.
558 : *
559 : * @return a dataset. If not NULL, it must be closed with GDALClose().
560 : *
561 : * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used
562 : * by gdalwarp initially, but is no longer needed.
563 : */
564 0 : GDALDatasetH GDALOpenVerticalShiftGrid(const char *pszProj4Geoidgrids,
565 : int *pbError)
566 : {
567 0 : char **papszGrids = CSLTokenizeString2(pszProj4Geoidgrids, ",", 0);
568 0 : const int nGridCount = CSLCount(papszGrids);
569 0 : if (nGridCount == 1)
570 : {
571 0 : CSLDestroy(papszGrids);
572 :
573 0 : bool bMissingOk = false;
574 0 : if (*pszProj4Geoidgrids == '@')
575 : {
576 0 : pszProj4Geoidgrids++;
577 0 : bMissingOk = true;
578 : }
579 0 : const CPLString osFilename(GetProj4Filename(pszProj4Geoidgrids));
580 0 : const char *const papszOpenOptions[] = {
581 : "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES", nullptr};
582 : GDALDatasetH hDS =
583 0 : GDALOpenEx(osFilename, 0, nullptr, papszOpenOptions, nullptr);
584 0 : if (hDS == nullptr)
585 : {
586 0 : CPLDebug("GDAL", "Cannot find file corresponding to %s",
587 : pszProj4Geoidgrids);
588 : }
589 0 : if (pbError)
590 0 : *pbError = (!bMissingOk && hDS == nullptr);
591 0 : return hDS;
592 : }
593 :
594 0 : CPLStringList aosFilenames;
595 0 : for (int i = nGridCount - 1; i >= 0; i--)
596 : {
597 0 : const char *pszName = papszGrids[i];
598 0 : bool bMissingOk = false;
599 0 : if (*pszName == '@')
600 : {
601 0 : pszName++;
602 0 : bMissingOk = true;
603 : }
604 0 : const CPLString osFilename(GetProj4Filename(pszName));
605 : VSIStatBufL sStat;
606 0 : if (osFilename.empty() || VSIStatL(osFilename, &sStat) != 0)
607 : {
608 0 : CPLDebug("GDAL", "Cannot find file corresponding to %s", pszName);
609 0 : if (!bMissingOk)
610 : {
611 0 : if (pbError)
612 0 : *pbError = true;
613 0 : CSLDestroy(papszGrids);
614 0 : return nullptr;
615 : }
616 : }
617 : else
618 : {
619 0 : aosFilenames.AddString(osFilename);
620 : }
621 : }
622 :
623 0 : CSLDestroy(papszGrids);
624 :
625 0 : if (aosFilenames.empty())
626 : {
627 0 : if (pbError)
628 0 : *pbError = false;
629 0 : return nullptr;
630 : }
631 :
632 0 : char **papszArgv = nullptr;
633 0 : papszArgv = CSLAddString(papszArgv, "-resolution");
634 0 : papszArgv = CSLAddString(papszArgv, "highest");
635 0 : papszArgv = CSLAddString(papszArgv, "-vrtnodata");
636 0 : papszArgv = CSLAddString(papszArgv, "-inf");
637 0 : papszArgv = CSLAddString(papszArgv, "-oo");
638 : papszArgv =
639 0 : CSLAddString(papszArgv, "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES");
640 0 : GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(papszArgv, nullptr);
641 0 : CSLDestroy(papszArgv);
642 0 : GDALDatasetH hDS = GDALBuildVRT("", aosFilenames.size(), nullptr,
643 0 : aosFilenames.List(), psOptions, nullptr);
644 0 : GDALBuildVRTOptionsFree(psOptions);
645 0 : if (pbError)
646 0 : *pbError = hDS != nullptr;
647 0 : return hDS;
648 : }
|
