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(GDALGeoTransform >) const 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(GDALGeoTransform >) const
143 : {
144 2 : return m_poSrcDataset->GetGeoTransform(gt);
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 23 : GDALGeoTransform srcGT;
351 23 : if (GDALDataset::FromHandle(hSrcDataset)->GetGeoTransform(srcGT) != 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 20 : GDALGeoTransform gridGT;
390 20 : if (GDALDataset::FromHandle(hGridDataset)->GetGeoTransform(gridGT) !=
391 : CE_None)
392 : {
393 1 : CPLError(CE_Failure, CPLE_NotSupported,
394 : "Grid dataset has no geotransform.");
395 1 : return nullptr;
396 : }
397 :
398 19 : OGRSpatialReferenceH hGridSRS = GDALGetSpatialRef(hGridDataset);
399 19 : if (hGridSRS == nullptr)
400 : {
401 1 : CPLError(CE_Failure, CPLE_NotSupported,
402 : "Grid dataset has no projection.");
403 1 : return nullptr;
404 : }
405 18 : if (GDALGetRasterCount(hGridDataset) != 1)
406 : {
407 1 : CPLError(CE_Failure, CPLE_NotSupported,
408 : "Only single band grid dataset is supported.");
409 1 : return nullptr;
410 : }
411 :
412 17 : GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDataset, 1));
413 17 : const char *pszDataType = CSLFetchNameValue(papszOptions, "DATATYPE");
414 17 : if (pszDataType)
415 2 : eDT = GDALGetDataTypeByName(pszDataType);
416 17 : if (eDT == GDT_Unknown)
417 : {
418 1 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid DATATYPE=%s",
419 : pszDataType);
420 1 : return nullptr;
421 : }
422 :
423 16 : const int nSrcXSize = GDALGetRasterXSize(hSrcDataset);
424 16 : const int nSrcYSize = GDALGetRasterYSize(hSrcDataset);
425 :
426 16 : double dfWestLongitudeDeg = 0.0;
427 16 : double dfSouthLatitudeDeg = 0.0;
428 16 : double dfEastLongitudeDeg = 0.0;
429 16 : double dfNorthLatitudeDeg = 0.0;
430 16 : GDALComputeAreaOfInterest(&oSrcSRS, srcGT.data(), nSrcXSize, nSrcYSize,
431 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
432 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
433 :
434 32 : CPLStringList aosOptions;
435 16 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
436 1 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
437 : {
438 : aosOptions.SetNameValue(
439 : "AREA_OF_INTEREST",
440 : CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
441 : dfSouthLatitudeDeg, dfEastLongitudeDeg,
442 15 : dfNorthLatitudeDeg));
443 : }
444 32 : void *hTransform = GDALCreateGenImgProjTransformer4(
445 16 : hGridSRS, gridGT.data(), OGRSpatialReference::ToHandle(&oSrcSRS),
446 16 : srcGT.data(), aosOptions.List());
447 16 : if (hTransform == nullptr)
448 3 : return nullptr;
449 13 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
450 13 : psWO->hSrcDS = hGridDataset;
451 13 : psWO->eResampleAlg = GRA_Bilinear;
452 13 : const char *pszResampling = CSLFetchNameValue(papszOptions, "RESAMPLING");
453 13 : if (pszResampling)
454 : {
455 3 : if (EQUAL(pszResampling, "NEAREST"))
456 1 : psWO->eResampleAlg = GRA_NearestNeighbour;
457 2 : else if (EQUAL(pszResampling, "BILINEAR"))
458 1 : psWO->eResampleAlg = GRA_Bilinear;
459 1 : else if (EQUAL(pszResampling, "CUBIC"))
460 1 : psWO->eResampleAlg = GRA_Cubic;
461 : }
462 13 : psWO->eWorkingDataType = GDT_Float32;
463 13 : int bHasNoData = FALSE;
464 13 : const double dfSrcNoData = GDALGetRasterNoDataValue(
465 : GDALGetRasterBand(hGridDataset, 1), &bHasNoData);
466 13 : if (bHasNoData)
467 : {
468 2 : psWO->padfSrcNoDataReal =
469 2 : static_cast<double *>(CPLMalloc(sizeof(double)));
470 2 : psWO->padfSrcNoDataReal[0] = dfSrcNoData;
471 : }
472 :
473 13 : psWO->padfDstNoDataReal = static_cast<double *>(CPLMalloc(sizeof(double)));
474 : const bool bErrorOnMissingShift =
475 13 : CPLFetchBool(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false);
476 13 : psWO->padfDstNoDataReal[0] =
477 13 : (bErrorOnMissingShift) ? -std::numeric_limits<float>::infinity() : 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), dfSrcUnitToMeter, dfDstUnitToMeter,
511 : // Undocumented option. For testing only
512 13 : atoi(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", "256")));
513 :
514 13 : poReprojectedGrid->ReleaseRef();
515 :
516 13 : if (!poOutDS->IsInitOK())
517 : {
518 1 : delete poOutDS;
519 1 : return nullptr;
520 : }
521 12 : poOutDS->SetDescription(GDALGetDescription(hSrcDataset));
522 12 : return reinterpret_cast<GDALDatasetH>(poOutDS);
523 : }
524 :
525 : /************************************************************************/
526 : /* GetProj4Filename() */
527 : /************************************************************************/
528 :
529 0 : static CPLString GetProj4Filename(const char *pszFilename)
530 : {
531 0 : CPLString osFilename;
532 :
533 : /* or fixed path: /name, ./name or ../name */
534 0 : if (!CPLIsFilenameRelative(pszFilename) || *pszFilename == '.')
535 : {
536 0 : return pszFilename;
537 : }
538 :
539 0 : PJ_GRID_INFO info = proj_grid_info(pszFilename);
540 0 : if (info.filename[0])
541 : {
542 0 : osFilename = info.filename;
543 : }
544 :
545 0 : return osFilename;
546 : }
547 :
548 : /************************************************************************/
549 : /* GDALOpenVerticalShiftGrid() */
550 : /************************************************************************/
551 :
552 : /** Load proj.4 geoidgrids as GDAL dataset
553 : *
554 : * @param pszProj4Geoidgrids Value of proj.4 geoidgrids parameter.
555 : * @param pbError If not NULL, the pointed value will be set to TRUE if an
556 : * error occurred.
557 : *
558 : * @return a dataset. If not NULL, it must be closed with GDALClose().
559 : *
560 : * @since GDAL 2.2
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 : }
|
