• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

OSGeo / gdal / 15885686134

25 Jun 2025 07:44PM UTC coverage: 71.084%. Remained the same
15885686134

push

github

rouault
gdal_priv.h: fix C++11 compatibility

573814 of 807237 relevant lines covered (71.08%)

250621.56 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

74.8
/alg/gdalapplyverticalshiftgrid.cpp
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 &gt) 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
GDALApplyVSGDataset::GDALApplyVSGDataset(GDALDataset *poSrcDataset,
13✔
86
                                         GDALDataset *poReprojectedGrid,
87
                                         GDALDataType eDT, bool bInverse,
88
                                         double dfSrcUnitToMeter,
89
                                         double dfDstUnitToMeter,
90
                                         int nBlockSize)
13✔
91
    : m_poSrcDataset(poSrcDataset), m_poReprojectedGrid(poReprojectedGrid),
92
      m_bInverse(bInverse), m_dfSrcUnitToMeter(dfSrcUnitToMeter),
93
      m_dfDstUnitToMeter(dfDstUnitToMeter)
13✔
94
{
95
    m_poSrcDataset->Reference();
13✔
96
    m_poReprojectedGrid->Reference();
13✔
97

98
    nRasterXSize = poSrcDataset->GetRasterXSize();
13✔
99
    nRasterYSize = poSrcDataset->GetRasterYSize();
13✔
100
    SetBand(1, new GDALApplyVSGRasterBand(eDT, nBlockSize));
13✔
101
}
13✔
102

103
/************************************************************************/
104
/*                       ~GDALApplyVSGDataset()                         */
105
/************************************************************************/
106

107
GDALApplyVSGDataset::~GDALApplyVSGDataset()
26✔
108
{
109
    GDALApplyVSGDataset::CloseDependentDatasets();
13✔
110
}
26✔
111

112
/************************************************************************/
113
/*                     CloseDependentDatasets()                         */
114
/************************************************************************/
115

116
int GDALApplyVSGDataset::CloseDependentDatasets()
13✔
117
{
118
    bool bRet = false;
13✔
119
    if (m_poSrcDataset != nullptr)
13✔
120
    {
121
        if (m_poSrcDataset->ReleaseRef())
13✔
122
        {
123
            bRet = true;
8✔
124
        }
125
        m_poSrcDataset = nullptr;
13✔
126
    }
127
    if (m_poReprojectedGrid != nullptr)
13✔
128
    {
129
        if (m_poReprojectedGrid->ReleaseRef())
13✔
130
        {
131
            bRet = true;
13✔
132
        }
133
        m_poReprojectedGrid = nullptr;
13✔
134
    }
135
    return bRet;
13✔
136
}
137

138
/************************************************************************/
139
/*                          GetGeoTransform()                           */
140
/************************************************************************/
141

142
CPLErr GDALApplyVSGDataset::GetGeoTransform(GDALGeoTransform &gt) const
2✔
143
{
144
    return m_poSrcDataset->GetGeoTransform(gt);
2✔
145
}
146

147
/************************************************************************/
148
/*                          GetSpatialRef()                             */
149
/************************************************************************/
150

151
const OGRSpatialReference *GDALApplyVSGDataset::GetSpatialRef() const
2✔
152
{
153
    return m_poSrcDataset->GetSpatialRef();
2✔
154
}
155

156
/************************************************************************/
157
/*                             IsInitOK()                               */
158
/************************************************************************/
159

160
bool GDALApplyVSGDataset::IsInitOK()
13✔
161
{
162
    GDALApplyVSGRasterBand *poBand =
163
        reinterpret_cast<GDALApplyVSGRasterBand *>(GetRasterBand(1));
13✔
164
    return poBand->m_pafSrcData != nullptr && poBand->m_pafGridData != nullptr;
13✔
165
}
166

167
/************************************************************************/
168
/*                       GDALApplyVSGRasterBand()                       */
169
/************************************************************************/
170

171
GDALApplyVSGRasterBand::GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize)
13✔
172
{
173
    eDataType = eDT;
13✔
174
    nBlockXSize = nBlockSize;
13✔
175
    nBlockYSize = nBlockSize;
13✔
176
    m_pafSrcData = static_cast<float *>(
13✔
177
        VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));
13✔
178
    m_pafGridData = static_cast<float *>(
13✔
179
        VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));
13✔
180
}
13✔
181

182
/************************************************************************/
183
/*                      ~GDALApplyVSGRasterBand()                       */
184
/************************************************************************/
185

186
GDALApplyVSGRasterBand::~GDALApplyVSGRasterBand()
26✔
187
{
188
    VSIFree(m_pafSrcData);
13✔
189
    VSIFree(m_pafGridData);
13✔
190
}
26✔
191

192
/************************************************************************/
193
/*                           GetNoDataValue()                           */
194
/************************************************************************/
195

196
double GDALApplyVSGRasterBand::GetNoDataValue(int *pbSuccess)
19✔
197
{
198
    GDALApplyVSGDataset *poGDS = reinterpret_cast<GDALApplyVSGDataset *>(poDS);
19✔
199
    return poGDS->m_poSrcDataset->GetRasterBand(1)->GetNoDataValue(pbSuccess);
19✔
200
}
201

202
/************************************************************************/
203
/*                              IReadBlock()                            */
204
/************************************************************************/
205

206
CPLErr GDALApplyVSGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
17✔
207
                                          void *pData)
208
{
209
    GDALApplyVSGDataset *poGDS = reinterpret_cast<GDALApplyVSGDataset *>(poDS);
17✔
210

211
    const int nXOff = nBlockXOff * nBlockXSize;
17✔
212
    const int nReqXSize = (nXOff > nRasterXSize - nBlockXSize)
34✔
213
                              ? nRasterXSize - nXOff
17✔
214
                              : nBlockXSize;
215
    const int nYOff = nBlockYOff * nBlockYSize;
17✔
216
    const int nReqYSize = (nYOff > nRasterYSize - nBlockYSize)
34✔
217
                              ? nRasterYSize - nYOff
17✔
218
                              : nBlockYSize;
219

220
    CPLErr eErr = poGDS->m_poSrcDataset->GetRasterBand(1)->RasterIO(
17✔
221
        GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafSrcData, nReqXSize,
17✔
222
        nReqYSize, GDT_Float32, sizeof(float), nBlockXSize * sizeof(float),
17✔
223
        nullptr);
224
    if (eErr == CE_None)
17✔
225
        eErr = poGDS->m_poReprojectedGrid->GetRasterBand(1)->RasterIO(
34✔
226
            GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafGridData,
17✔
227
            nReqXSize, nReqYSize, GDT_Float32, sizeof(float),
228
            nBlockXSize * sizeof(float), nullptr);
17✔
229
    if (eErr == CE_None)
17✔
230
    {
231
        const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
17✔
232
        int bHasNoData = FALSE;
17✔
233
        float fNoDataValue = static_cast<float>(GetNoDataValue(&bHasNoData));
17✔
234
        for (int iY = 0; iY < nReqYSize; iY++)
279✔
235
        {
236
            for (int iX = 0; iX < nReqXSize; iX++)
4,666✔
237
            {
238
                const float fSrcVal = m_pafSrcData[iY * nBlockXSize + iX];
4,404✔
239
                const float fGridVal = m_pafGridData[iY * nBlockXSize + iX];
4,404✔
240
                if (bHasNoData && fSrcVal == fNoDataValue)
4,404✔
241
                {
242
                }
243
                else if (std::isinf(fGridVal))
4,403✔
244
                {
245
                    CPLError(CE_Failure, CPLE_AppDefined,
2✔
246
                             "Missing vertical grid value at source (%d,%d)",
247
                             nXOff + iX, nYOff + iY);
248
                    return CE_Failure;
2✔
249
                }
250
                else if (poGDS->m_bInverse)
4,401✔
251
                {
252
                    m_pafSrcData[iY * nBlockXSize + iX] = static_cast<float>(
800✔
253
                        (fSrcVal * poGDS->m_dfSrcUnitToMeter - fGridVal) /
800✔
254
                        poGDS->m_dfDstUnitToMeter);
800✔
255
                }
256
                else
257
                {
258
                    m_pafSrcData[iY * nBlockXSize + iX] = static_cast<float>(
3,601✔
259
                        (fSrcVal * poGDS->m_dfSrcUnitToMeter + fGridVal) /
3,601✔
260
                        poGDS->m_dfDstUnitToMeter);
3,601✔
261
                }
262
            }
263
            GDALCopyWords(
262✔
264
                m_pafSrcData + iY * nBlockXSize, GDT_Float32, sizeof(float),
262✔
265
                static_cast<GByte *>(pData) + iY * nBlockXSize * nDTSize,
262✔
266
                eDataType, nDTSize, nReqXSize);
267
        }
268
    }
269

270
    return eErr;
15✔
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
GDALDatasetH GDALApplyVerticalShiftGrid(GDALDatasetH hSrcDataset,
23✔
342
                                        GDALDatasetH hGridDataset, int bInverse,
343
                                        double dfSrcUnitToMeter,
344
                                        double dfDstUnitToMeter,
345
                                        const char *const *papszOptions)
346
{
347
    VALIDATE_POINTER1(hSrcDataset, "GDALApplyVerticalShiftGrid", nullptr);
23✔
348
    VALIDATE_POINTER1(hGridDataset, "GDALApplyVerticalShiftGrid", nullptr);
23✔
349

350
    GDALGeoTransform srcGT;
23✔
351
    if (GDALDataset::FromHandle(hSrcDataset)->GetGeoTransform(srcGT) != CE_None)
23✔
352
    {
353
        CPLError(CE_Failure, CPLE_NotSupported,
1✔
354
                 "Source dataset has no geotransform.");
355
        return nullptr;
1✔
356
    }
357
    const char *pszSrcProjection = CSLFetchNameValue(papszOptions, "SRC_SRS");
22✔
358
    OGRSpatialReference oSrcSRS;
44✔
359
    if (pszSrcProjection != nullptr && pszSrcProjection[0] != '\0')
22✔
360
    {
361
        oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
×
362
        oSrcSRS.SetFromUserInput(pszSrcProjection);
×
363
    }
364
    else
365
    {
366
        auto poSRS = GDALDataset::FromHandle(hSrcDataset)->GetSpatialRef();
22✔
367
        if (poSRS)
22✔
368
            oSrcSRS = *poSRS;
21✔
369
    }
370

371
    if (oSrcSRS.IsCompound())
22✔
372
    {
373
        oSrcSRS.StripVertical();
×
374
    }
375

376
    if (oSrcSRS.IsEmpty())
22✔
377
    {
378
        CPLError(CE_Failure, CPLE_NotSupported,
1✔
379
                 "Source dataset has no projection.");
380
        return nullptr;
1✔
381
    }
382
    if (GDALGetRasterCount(hSrcDataset) != 1)
21✔
383
    {
384
        CPLError(CE_Failure, CPLE_NotSupported,
1✔
385
                 "Only single band source dataset is supported.");
386
        return nullptr;
1✔
387
    }
388

389
    GDALGeoTransform gridGT;
20✔
390
    if (GDALDataset::FromHandle(hGridDataset)->GetGeoTransform(gridGT) !=
20✔
391
        CE_None)
392
    {
393
        CPLError(CE_Failure, CPLE_NotSupported,
1✔
394
                 "Grid dataset has no geotransform.");
395
        return nullptr;
1✔
396
    }
397

398
    OGRSpatialReferenceH hGridSRS = GDALGetSpatialRef(hGridDataset);
19✔
399
    if (hGridSRS == nullptr)
19✔
400
    {
401
        CPLError(CE_Failure, CPLE_NotSupported,
1✔
402
                 "Grid dataset has no projection.");
403
        return nullptr;
1✔
404
    }
405
    if (GDALGetRasterCount(hGridDataset) != 1)
18✔
406
    {
407
        CPLError(CE_Failure, CPLE_NotSupported,
1✔
408
                 "Only single band grid dataset is supported.");
409
        return nullptr;
1✔
410
    }
411

412
    GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDataset, 1));
17✔
413
    const char *pszDataType = CSLFetchNameValue(papszOptions, "DATATYPE");
17✔
414
    if (pszDataType)
17✔
415
        eDT = GDALGetDataTypeByName(pszDataType);
2✔
416
    if (eDT == GDT_Unknown)
17✔
417
    {
418
        CPLError(CE_Failure, CPLE_NotSupported, "Invalid DATATYPE=%s",
1✔
419
                 pszDataType);
420
        return nullptr;
1✔
421
    }
422

423
    const int nSrcXSize = GDALGetRasterXSize(hSrcDataset);
16✔
424
    const int nSrcYSize = GDALGetRasterYSize(hSrcDataset);
16✔
425

426
    double dfWestLongitudeDeg = 0.0;
16✔
427
    double dfSouthLatitudeDeg = 0.0;
16✔
428
    double dfEastLongitudeDeg = 0.0;
16✔
429
    double dfNorthLatitudeDeg = 0.0;
16✔
430
    GDALComputeAreaOfInterest(&oSrcSRS, srcGT.data(), nSrcXSize, nSrcYSize,
16✔
431
                              dfWestLongitudeDeg, dfSouthLatitudeDeg,
432
                              dfEastLongitudeDeg, dfNorthLatitudeDeg);
433

434
    CPLStringList aosOptions;
32✔
435
    if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
16✔
436
          dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
1✔
437
    {
438
        aosOptions.SetNameValue(
439
            "AREA_OF_INTEREST",
440
            CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
441
                       dfSouthLatitudeDeg, dfEastLongitudeDeg,
442
                       dfNorthLatitudeDeg));
15✔
443
    }
444
    void *hTransform = GDALCreateGenImgProjTransformer4(
32✔
445
        hGridSRS, gridGT.data(), OGRSpatialReference::ToHandle(&oSrcSRS),
16✔
446
        srcGT.data(), aosOptions.List());
16✔
447
    if (hTransform == nullptr)
16✔
448
        return nullptr;
3✔
449
    GDALWarpOptions *psWO = GDALCreateWarpOptions();
13✔
450
    psWO->hSrcDS = hGridDataset;
13✔
451
    psWO->eResampleAlg = GRA_Bilinear;
13✔
452
    const char *pszResampling = CSLFetchNameValue(papszOptions, "RESAMPLING");
13✔
453
    if (pszResampling)
13✔
454
    {
455
        if (EQUAL(pszResampling, "NEAREST"))
3✔
456
            psWO->eResampleAlg = GRA_NearestNeighbour;
1✔
457
        else if (EQUAL(pszResampling, "BILINEAR"))
2✔
458
            psWO->eResampleAlg = GRA_Bilinear;
1✔
459
        else if (EQUAL(pszResampling, "CUBIC"))
1✔
460
            psWO->eResampleAlg = GRA_Cubic;
1✔
461
    }
462
    psWO->eWorkingDataType = GDT_Float32;
13✔
463
    int bHasNoData = FALSE;
13✔
464
    const double dfSrcNoData = GDALGetRasterNoDataValue(
13✔
465
        GDALGetRasterBand(hGridDataset, 1), &bHasNoData);
466
    if (bHasNoData)
13✔
467
    {
468
        psWO->padfSrcNoDataReal =
2✔
469
            static_cast<double *>(CPLMalloc(sizeof(double)));
2✔
470
        psWO->padfSrcNoDataReal[0] = dfSrcNoData;
2✔
471
    }
472

473
    psWO->padfDstNoDataReal = static_cast<double *>(CPLMalloc(sizeof(double)));
13✔
474
    const bool bErrorOnMissingShift =
475
        CPLFetchBool(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false);
13✔
476
    psWO->padfDstNoDataReal[0] =
13✔
477
        (bErrorOnMissingShift) ? -std::numeric_limits<float>::infinity() : 0.0;
13✔
478
    psWO->papszWarpOptions =
13✔
479
        CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
13✔
480

481
    psWO->pfnTransformer = GDALGenImgProjTransform;
13✔
482
    psWO->pTransformerArg = hTransform;
13✔
483
    const double dfMaxError =
484
        CPLAtof(CSLFetchNameValueDef(papszOptions, "MAX_ERROR", "0.125"));
13✔
485
    if (dfMaxError > 0.0)
13✔
486
    {
487
        psWO->pTransformerArg = GDALCreateApproxTransformer(
13✔
488
            psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);
489
        psWO->pfnTransformer = GDALApproxTransform;
13✔
490
        GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);
13✔
491
    }
492
    psWO->nBandCount = 1;
13✔
493
    psWO->panSrcBands = static_cast<int *>(CPLMalloc(sizeof(int)));
13✔
494
    psWO->panSrcBands[0] = 1;
13✔
495
    psWO->panDstBands = static_cast<int *>(CPLMalloc(sizeof(int)));
13✔
496
    psWO->panDstBands[0] = 1;
13✔
497

498
    VRTWarpedDataset *poReprojectedGrid =
499
        new VRTWarpedDataset(nSrcXSize, nSrcYSize);
13✔
500
    // This takes a reference on hGridDataset
501
    CPLErr eErr = poReprojectedGrid->Initialize(psWO);
13✔
502
    CPLAssert(eErr == CE_None);
13✔
503
    CPL_IGNORE_RET_VAL(eErr);
13✔
504
    GDALDestroyWarpOptions(psWO);
13✔
505
    poReprojectedGrid->SetGeoTransform(srcGT);
13✔
506
    poReprojectedGrid->AddBand(GDT_Float32, nullptr);
13✔
507

508
    GDALApplyVSGDataset *poOutDS = new GDALApplyVSGDataset(
509
        GDALDataset::FromHandle(hSrcDataset), poReprojectedGrid, eDT,
13✔
510
        CPL_TO_BOOL(bInverse), dfSrcUnitToMeter, dfDstUnitToMeter,
13✔
511
        // Undocumented option. For testing only
512
        atoi(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", "256")));
13✔
513

514
    poReprojectedGrid->ReleaseRef();
13✔
515

516
    if (!poOutDS->IsInitOK())
13✔
517
    {
518
        delete poOutDS;
1✔
519
        return nullptr;
1✔
520
    }
521
    poOutDS->SetDescription(GDALGetDescription(hSrcDataset));
12✔
522
    return reinterpret_cast<GDALDatasetH>(poOutDS);
12✔
523
}
524

525
/************************************************************************/
526
/*                           GetProj4Filename()                         */
527
/************************************************************************/
528

529
static CPLString GetProj4Filename(const char *pszFilename)
×
530
{
531
    CPLString osFilename;
×
532

533
    /* or fixed path: /name, ./name or ../name */
534
    if (!CPLIsFilenameRelative(pszFilename) || *pszFilename == '.')
×
535
    {
536
        return pszFilename;
×
537
    }
538

539
    PJ_GRID_INFO info = proj_grid_info(pszFilename);
×
540
    if (info.filename[0])
×
541
    {
542
        osFilename = info.filename;
×
543
    }
544

545
    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
GDALDatasetH GDALOpenVerticalShiftGrid(const char *pszProj4Geoidgrids,
×
565
                                       int *pbError)
566
{
567
    char **papszGrids = CSLTokenizeString2(pszProj4Geoidgrids, ",", 0);
×
568
    const int nGridCount = CSLCount(papszGrids);
×
569
    if (nGridCount == 1)
×
570
    {
571
        CSLDestroy(papszGrids);
×
572

573
        bool bMissingOk = false;
×
574
        if (*pszProj4Geoidgrids == '@')
×
575
        {
576
            pszProj4Geoidgrids++;
×
577
            bMissingOk = true;
×
578
        }
579
        const CPLString osFilename(GetProj4Filename(pszProj4Geoidgrids));
×
580
        const char *const papszOpenOptions[] = {
×
581
            "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES", nullptr};
582
        GDALDatasetH hDS =
583
            GDALOpenEx(osFilename, 0, nullptr, papszOpenOptions, nullptr);
×
584
        if (hDS == nullptr)
×
585
        {
586
            CPLDebug("GDAL", "Cannot find file corresponding to %s",
×
587
                     pszProj4Geoidgrids);
588
        }
589
        if (pbError)
×
590
            *pbError = (!bMissingOk && hDS == nullptr);
×
591
        return hDS;
×
592
    }
593

594
    CPLStringList aosFilenames;
×
595
    for (int i = nGridCount - 1; i >= 0; i--)
×
596
    {
597
        const char *pszName = papszGrids[i];
×
598
        bool bMissingOk = false;
×
599
        if (*pszName == '@')
×
600
        {
601
            pszName++;
×
602
            bMissingOk = true;
×
603
        }
604
        const CPLString osFilename(GetProj4Filename(pszName));
×
605
        VSIStatBufL sStat;
606
        if (osFilename.empty() || VSIStatL(osFilename, &sStat) != 0)
×
607
        {
608
            CPLDebug("GDAL", "Cannot find file corresponding to %s", pszName);
×
609
            if (!bMissingOk)
×
610
            {
611
                if (pbError)
×
612
                    *pbError = true;
×
613
                CSLDestroy(papszGrids);
×
614
                return nullptr;
×
615
            }
616
        }
617
        else
618
        {
619
            aosFilenames.AddString(osFilename);
×
620
        }
621
    }
622

623
    CSLDestroy(papszGrids);
×
624

625
    if (aosFilenames.empty())
×
626
    {
627
        if (pbError)
×
628
            *pbError = false;
×
629
        return nullptr;
×
630
    }
631

632
    char **papszArgv = nullptr;
×
633
    papszArgv = CSLAddString(papszArgv, "-resolution");
×
634
    papszArgv = CSLAddString(papszArgv, "highest");
×
635
    papszArgv = CSLAddString(papszArgv, "-vrtnodata");
×
636
    papszArgv = CSLAddString(papszArgv, "-inf");
×
637
    papszArgv = CSLAddString(papszArgv, "-oo");
×
638
    papszArgv =
639
        CSLAddString(papszArgv, "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES");
×
640
    GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(papszArgv, nullptr);
×
641
    CSLDestroy(papszArgv);
×
642
    GDALDatasetH hDS = GDALBuildVRT("", aosFilenames.size(), nullptr,
×
643
                                    aosFilenames.List(), psOptions, nullptr);
×
644
    GDALBuildVRTOptionsFree(psOptions);
×
645
    if (pbError)
×
646
        *pbError = hDS != nullptr;
×
647
    return hDS;
×
648
}
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc