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