Coverage Report

Created: 2025-06-13 06:18

/src/gdal/ogr/ogrsf_frmts/generic/ogrwarpedlayer.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  OpenGIS Simple Features Reference Implementation
4
 * Purpose:  Implements OGRWarpedLayer class
5
 * Author:   Even Rouault, even dot rouault at spatialys.com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2012-2014, Even Rouault <even dot rouault at spatialys.com>
9
 *
10
 * SPDX-License-Identifier: MIT
11
 ****************************************************************************/
12
13
#ifndef DOXYGEN_SKIP
14
15
#include <cmath>
16
17
#include "ogrwarpedlayer.h"
18
19
/************************************************************************/
20
/*                          OGRWarpedLayer()                            */
21
/************************************************************************/
22
23
OGRWarpedLayer::OGRWarpedLayer(OGRLayer *poDecoratedLayer, int iGeomField,
24
                               int bTakeOwnership,
25
                               OGRCoordinateTransformation *poCT,
26
                               OGRCoordinateTransformation *poReversedCT)
27
0
    : OGRLayerDecorator(poDecoratedLayer, bTakeOwnership),
28
0
      m_poFeatureDefn(nullptr), m_iGeomField(iGeomField), m_poCT(poCT),
29
0
      m_poReversedCT(poReversedCT),
30
0
      m_poSRS(const_cast<OGRSpatialReference *>(m_poCT->GetTargetCS()))
31
0
{
32
0
    CPLAssert(poCT != nullptr);
33
0
    SetDescription(poDecoratedLayer->GetDescription());
34
35
0
    if (m_poSRS != nullptr)
36
0
    {
37
0
        m_poSRS->Reference();
38
0
    }
39
0
}
Unexecuted instantiation: OGRWarpedLayer::OGRWarpedLayer(OGRLayer*, int, int, OGRCoordinateTransformation*, OGRCoordinateTransformation*)
Unexecuted instantiation: OGRWarpedLayer::OGRWarpedLayer(OGRLayer*, int, int, OGRCoordinateTransformation*, OGRCoordinateTransformation*)
40
41
/************************************************************************/
42
/*                         ~OGRWarpedLayer()                            */
43
/************************************************************************/
44
45
OGRWarpedLayer::~OGRWarpedLayer()
46
0
{
47
0
    if (m_poFeatureDefn != nullptr)
48
0
        m_poFeatureDefn->Release();
49
0
    if (m_poSRS != nullptr)
50
0
        m_poSRS->Release();
51
0
    delete m_poCT;
52
0
    delete m_poReversedCT;
53
0
}
54
55
/************************************************************************/
56
/*                         ISetSpatialFilter()                          */
57
/************************************************************************/
58
59
OGRErr OGRWarpedLayer::ISetSpatialFilter(int iGeomField,
60
                                         const OGRGeometry *poGeom)
61
0
{
62
63
0
    m_iGeomFieldFilter = iGeomField;
64
0
    if (InstallFilter(poGeom))
65
0
        ResetReading();
66
67
0
    if (m_iGeomFieldFilter == m_iGeomField)
68
0
    {
69
0
        if (poGeom == nullptr || m_poReversedCT == nullptr)
70
0
        {
71
0
            return m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter,
72
0
                                                        nullptr);
73
0
        }
74
0
        else
75
0
        {
76
0
            OGREnvelope sEnvelope;
77
0
            poGeom->getEnvelope(&sEnvelope);
78
0
            if (std::isinf(sEnvelope.MinX) && std::isinf(sEnvelope.MinY) &&
79
0
                std::isinf(sEnvelope.MaxX) && std::isinf(sEnvelope.MaxY))
80
0
            {
81
0
                return m_poDecoratedLayer->SetSpatialFilterRect(
82
0
                    m_iGeomFieldFilter, sEnvelope.MinX, sEnvelope.MinY,
83
0
                    sEnvelope.MaxX, sEnvelope.MaxY);
84
0
            }
85
0
            else if (ReprojectEnvelope(&sEnvelope, m_poReversedCT))
86
0
            {
87
0
                return m_poDecoratedLayer->SetSpatialFilterRect(
88
0
                    m_iGeomFieldFilter, sEnvelope.MinX, sEnvelope.MinY,
89
0
                    sEnvelope.MaxX, sEnvelope.MaxY);
90
0
            }
91
0
            else
92
0
            {
93
0
                return m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter,
94
0
                                                            nullptr);
95
0
            }
96
0
        }
97
0
    }
98
0
    else
99
0
    {
100
0
        return m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter, poGeom);
101
0
    }
102
0
}
103
104
/************************************************************************/
105
/*                         TranslateFeature()                           */
106
/************************************************************************/
107
108
void OGRWarpedLayer::TranslateFeature(
109
    std::unique_ptr<OGRFeature> poSrcFeature,
110
    std::vector<std::unique_ptr<OGRFeature>> &apoOutFeatures)
111
0
{
112
0
    apoOutFeatures.push_back(
113
0
        SrcFeatureToWarpedFeature(std::move(poSrcFeature)));
114
0
}
115
116
/************************************************************************/
117
/*                     SrcFeatureToWarpedFeature()                      */
118
/************************************************************************/
119
120
std::unique_ptr<OGRFeature>
121
OGRWarpedLayer::SrcFeatureToWarpedFeature(std::unique_ptr<OGRFeature> poFeature)
122
0
{
123
    // This is safe to do here as they have matching attribute and geometry
124
    // fields
125
0
    poFeature->SetFDefnUnsafe(GetLayerDefn());
126
127
0
    OGRGeometry *poGeom = poFeature->GetGeomFieldRef(m_iGeomField);
128
0
    if (poGeom && poGeom->transform(m_poCT) != OGRERR_NONE)
129
0
    {
130
0
        delete poFeature->StealGeometry(m_iGeomField);
131
0
    }
132
133
0
    return poFeature;
134
0
}
135
136
/************************************************************************/
137
/*                     WarpedFeatureToSrcFeature()                      */
138
/************************************************************************/
139
140
std::unique_ptr<OGRFeature>
141
OGRWarpedLayer::WarpedFeatureToSrcFeature(std::unique_ptr<OGRFeature> poFeature)
142
0
{
143
    // This is safe to do here as they have matching attribute and geometry
144
    // fields
145
0
    poFeature->SetFDefnUnsafe(m_poDecoratedLayer->GetLayerDefn());
146
147
0
    OGRGeometry *poGeom = poFeature->GetGeomFieldRef(m_iGeomField);
148
0
    if (poGeom &&
149
0
        (!m_poReversedCT || poGeom->transform(m_poReversedCT) != OGRERR_NONE))
150
0
    {
151
0
        return nullptr;
152
0
    }
153
154
0
    return poFeature;
155
0
}
156
157
/************************************************************************/
158
/*                          GetNextFeature()                            */
159
/************************************************************************/
160
161
OGRFeature *OGRWarpedLayer::GetNextFeature()
162
0
{
163
0
    while (true)
164
0
    {
165
0
        auto poFeature =
166
0
            std::unique_ptr<OGRFeature>(m_poDecoratedLayer->GetNextFeature());
167
0
        if (!poFeature)
168
0
            return nullptr;
169
170
0
        auto poFeatureNew = SrcFeatureToWarpedFeature(std::move(poFeature));
171
0
        const OGRGeometry *poGeom = poFeatureNew->GetGeomFieldRef(m_iGeomField);
172
0
        if (m_poFilterGeom != nullptr && !FilterGeometry(poGeom))
173
0
        {
174
0
            continue;
175
0
        }
176
177
0
        return poFeatureNew.release();
178
0
    }
179
0
}
180
181
/************************************************************************/
182
/*                             GetFeature()                             */
183
/************************************************************************/
184
185
OGRFeature *OGRWarpedLayer::GetFeature(GIntBig nFID)
186
0
{
187
0
    auto poFeature =
188
0
        std::unique_ptr<OGRFeature>(m_poDecoratedLayer->GetFeature(nFID));
189
0
    if (poFeature)
190
0
    {
191
0
        poFeature = SrcFeatureToWarpedFeature(std::move(poFeature));
192
0
    }
193
0
    return poFeature.release();
194
0
}
195
196
/************************************************************************/
197
/*                             ISetFeature()                             */
198
/************************************************************************/
199
200
OGRErr OGRWarpedLayer::ISetFeature(OGRFeature *poFeature)
201
0
{
202
0
    auto poFeatureNew = WarpedFeatureToSrcFeature(
203
0
        std::unique_ptr<OGRFeature>(poFeature->Clone()));
204
0
    if (!poFeatureNew)
205
0
        return OGRERR_FAILURE;
206
207
0
    return m_poDecoratedLayer->SetFeature(poFeatureNew.get());
208
0
}
209
210
/************************************************************************/
211
/*                            ICreateFeature()                           */
212
/************************************************************************/
213
214
OGRErr OGRWarpedLayer::ICreateFeature(OGRFeature *poFeature)
215
0
{
216
0
    auto poFeatureNew = WarpedFeatureToSrcFeature(
217
0
        std::unique_ptr<OGRFeature>(poFeature->Clone()));
218
0
    if (!poFeatureNew)
219
0
        return OGRERR_FAILURE;
220
221
0
    return m_poDecoratedLayer->CreateFeature(poFeatureNew.get());
222
0
}
223
224
/************************************************************************/
225
/*                            IUpsertFeature()                          */
226
/************************************************************************/
227
228
OGRErr OGRWarpedLayer::IUpsertFeature(OGRFeature *poFeature)
229
0
{
230
0
    auto poFeatureNew = WarpedFeatureToSrcFeature(
231
0
        std::unique_ptr<OGRFeature>(poFeature->Clone()));
232
0
    if (poFeatureNew == nullptr)
233
0
        return OGRERR_FAILURE;
234
235
0
    return m_poDecoratedLayer->UpsertFeature(poFeatureNew.get());
236
0
}
237
238
/************************************************************************/
239
/*                            IUpdateFeature()                          */
240
/************************************************************************/
241
242
OGRErr OGRWarpedLayer::IUpdateFeature(OGRFeature *poFeature,
243
                                      int nUpdatedFieldsCount,
244
                                      const int *panUpdatedFieldsIdx,
245
                                      int nUpdatedGeomFieldsCount,
246
                                      const int *panUpdatedGeomFieldsIdx,
247
                                      bool bUpdateStyleString)
248
0
{
249
0
    auto poFeatureNew = WarpedFeatureToSrcFeature(
250
0
        std::unique_ptr<OGRFeature>(poFeature->Clone()));
251
0
    if (!poFeatureNew)
252
0
        return OGRERR_FAILURE;
253
254
0
    return m_poDecoratedLayer->UpdateFeature(
255
0
        poFeatureNew.get(), nUpdatedFieldsCount, panUpdatedFieldsIdx,
256
0
        nUpdatedGeomFieldsCount, panUpdatedGeomFieldsIdx, bUpdateStyleString);
257
0
}
258
259
/************************************************************************/
260
/*                            GetLayerDefn()                           */
261
/************************************************************************/
262
263
OGRFeatureDefn *OGRWarpedLayer::GetLayerDefn()
264
0
{
265
0
    if (m_poFeatureDefn != nullptr)
266
0
        return m_poFeatureDefn;
267
268
0
    m_poFeatureDefn = m_poDecoratedLayer->GetLayerDefn()->Clone();
269
0
    m_poFeatureDefn->Reference();
270
0
    if (m_poFeatureDefn->GetGeomFieldCount() > 0)
271
0
        m_poFeatureDefn->GetGeomFieldDefn(m_iGeomField)->SetSpatialRef(m_poSRS);
272
273
0
    return m_poFeatureDefn;
274
0
}
275
276
/************************************************************************/
277
/*                            GetSpatialRef()                           */
278
/************************************************************************/
279
280
OGRSpatialReference *OGRWarpedLayer::GetSpatialRef()
281
0
{
282
0
    if (m_iGeomField == 0)
283
0
        return m_poSRS;
284
0
    else
285
0
        return OGRLayer::GetSpatialRef();
286
0
}
287
288
/************************************************************************/
289
/*                           GetFeatureCount()                          */
290
/************************************************************************/
291
292
GIntBig OGRWarpedLayer::GetFeatureCount(int bForce)
293
0
{
294
0
    if (m_poFilterGeom == nullptr)
295
0
        return m_poDecoratedLayer->GetFeatureCount(bForce);
296
297
0
    return OGRLayer::GetFeatureCount(bForce);
298
0
}
299
300
/************************************************************************/
301
/*                             IGetExtent()                             */
302
/************************************************************************/
303
304
OGRErr OGRWarpedLayer::IGetExtent(int iGeomField, OGREnvelope *psExtent,
305
                                  bool bForce)
306
0
{
307
0
    if (iGeomField == m_iGeomField)
308
0
    {
309
0
        if (sStaticEnvelope.IsInit())
310
0
        {
311
0
            *psExtent = sStaticEnvelope;
312
0
            return OGRERR_NONE;
313
0
        }
314
315
0
        OGREnvelope sExtent;
316
0
        OGRErr eErr =
317
0
            m_poDecoratedLayer->GetExtent(m_iGeomField, &sExtent, bForce);
318
0
        if (eErr != OGRERR_NONE)
319
0
            return eErr;
320
321
0
        if (ReprojectEnvelope(&sExtent, m_poCT))
322
0
        {
323
0
            *psExtent = sExtent;
324
0
            return OGRERR_NONE;
325
0
        }
326
0
        else
327
0
            return OGRERR_FAILURE;
328
0
    }
329
0
    else
330
0
        return m_poDecoratedLayer->GetExtent(iGeomField, psExtent, bForce);
331
0
}
332
333
/************************************************************************/
334
/*                    TransformAndUpdateBBAndReturnX()                  */
335
/************************************************************************/
336
337
static double TransformAndUpdateBBAndReturnX(OGRCoordinateTransformation *poCT,
338
                                             double dfX, double dfY,
339
                                             double &dfMinX, double &dfMinY,
340
                                             double &dfMaxX, double &dfMaxY)
341
0
{
342
0
    int bSuccess = FALSE;
343
0
    poCT->Transform(1, &dfX, &dfY, nullptr, nullptr, &bSuccess);
344
0
    if (bSuccess)
345
0
    {
346
0
        if (dfX < dfMinX)
347
0
            dfMinX = dfX;
348
0
        if (dfY < dfMinY)
349
0
            dfMinY = dfY;
350
0
        if (dfX > dfMaxX)
351
0
            dfMaxX = dfX;
352
0
        if (dfY > dfMaxY)
353
0
            dfMaxY = dfY;
354
0
        return dfX;
355
0
    }
356
357
0
    return 0.0;
358
0
}
359
360
/************************************************************************/
361
/*                            FindXDiscontinuity()                      */
362
/************************************************************************/
363
364
static void FindXDiscontinuity(OGRCoordinateTransformation *poCT, double dfX1,
365
                               double dfX2, double dfY, double &dfMinX,
366
                               double &dfMinY, double &dfMaxX, double &dfMaxY,
367
                               int nRecLevel = 0)
368
0
{
369
0
    double dfXMid = (dfX1 + dfX2) / 2;
370
371
0
    double dfWrkX1 = TransformAndUpdateBBAndReturnX(poCT, dfX1, dfY, dfMinX,
372
0
                                                    dfMinY, dfMaxX, dfMaxY);
373
0
    double dfWrkXMid = TransformAndUpdateBBAndReturnX(poCT, dfXMid, dfY, dfMinX,
374
0
                                                      dfMinY, dfMaxX, dfMaxY);
375
0
    double dfWrkX2 = TransformAndUpdateBBAndReturnX(poCT, dfX2, dfY, dfMinX,
376
0
                                                    dfMinY, dfMaxX, dfMaxY);
377
378
0
    double dfDX1 = dfWrkXMid - dfWrkX1;
379
0
    double dfDX2 = dfWrkX2 - dfWrkXMid;
380
381
0
    if (dfDX1 * dfDX2 < 0 && nRecLevel < 30)
382
0
    {
383
0
        FindXDiscontinuity(poCT, dfX1, dfXMid, dfY, dfMinX, dfMinY, dfMaxX,
384
0
                           dfMaxY, nRecLevel + 1);
385
0
        FindXDiscontinuity(poCT, dfXMid, dfX2, dfY, dfMinX, dfMinY, dfMaxX,
386
0
                           dfMaxY, nRecLevel + 1);
387
0
    }
388
0
}
389
390
/************************************************************************/
391
/*                            ReprojectEnvelope()                       */
392
/************************************************************************/
393
394
int OGRWarpedLayer::ReprojectEnvelope(OGREnvelope *psEnvelope,
395
                                      OGRCoordinateTransformation *poCT)
396
0
{
397
0
    const int NSTEP = 20;
398
0
    double dfXStep = (psEnvelope->MaxX - psEnvelope->MinX) / NSTEP;
399
0
    double dfYStep = (psEnvelope->MaxY - psEnvelope->MinY) / NSTEP;
400
401
0
    double *padfX = static_cast<double *>(
402
0
        VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(double)));
403
0
    double *padfY = static_cast<double *>(
404
0
        VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(double)));
405
0
    int *pabSuccess = static_cast<int *>(
406
0
        VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(int)));
407
0
    if (padfX == nullptr || padfY == nullptr || pabSuccess == nullptr)
408
0
    {
409
0
        VSIFree(padfX);
410
0
        VSIFree(padfY);
411
0
        VSIFree(pabSuccess);
412
0
        return FALSE;
413
0
    }
414
415
0
    for (int j = 0; j <= NSTEP; j++)
416
0
    {
417
0
        for (int i = 0; i <= NSTEP; i++)
418
0
        {
419
0
            padfX[j * (NSTEP + 1) + i] = psEnvelope->MinX + i * dfXStep;
420
0
            padfY[j * (NSTEP + 1) + i] = psEnvelope->MinY + j * dfYStep;
421
0
        }
422
0
    }
423
424
0
    int bRet = FALSE;
425
426
0
    if (poCT->Transform((NSTEP + 1) * (NSTEP + 1), padfX, padfY, nullptr,
427
0
                        nullptr, pabSuccess))
428
0
    {
429
0
        double dfMinX = 0.0;
430
0
        double dfMinY = 0.0;
431
0
        double dfMaxX = 0.0;
432
0
        double dfMaxY = 0.0;
433
0
        int bSet = FALSE;
434
0
        for (int j = 0; j <= NSTEP; j++)
435
0
        {
436
0
            double dfXOld = 0.0;
437
0
            double dfDXOld = 0.0;
438
0
            int iOld = -1;
439
0
            int iOldOld = -1;
440
0
            for (int i = 0; i <= NSTEP; i++)
441
0
            {
442
0
                if (pabSuccess[j * (NSTEP + 1) + i])
443
0
                {
444
0
                    double dfX = padfX[j * (NSTEP + 1) + i];
445
0
                    double dfY = padfY[j * (NSTEP + 1) + i];
446
447
0
                    if (!bSet)
448
0
                    {
449
0
                        dfMinX = dfX;
450
0
                        dfMaxX = dfX;
451
0
                        dfMinY = dfY;
452
0
                        dfMaxY = dfY;
453
0
                        bSet = TRUE;
454
0
                    }
455
0
                    else
456
0
                    {
457
0
                        if (dfX < dfMinX)
458
0
                            dfMinX = dfX;
459
0
                        if (dfY < dfMinY)
460
0
                            dfMinY = dfY;
461
0
                        if (dfX > dfMaxX)
462
0
                            dfMaxX = dfX;
463
0
                        if (dfY > dfMaxY)
464
0
                            dfMaxY = dfY;
465
0
                    }
466
467
0
                    if (iOld >= 0)
468
0
                    {
469
0
                        double dfDXNew = dfX - dfXOld;
470
0
                        if (iOldOld >= 0 && dfDXNew * dfDXOld < 0)
471
0
                        {
472
0
                            FindXDiscontinuity(
473
0
                                poCT, psEnvelope->MinX + iOldOld * dfXStep,
474
0
                                psEnvelope->MinX + i * dfXStep,
475
0
                                psEnvelope->MinY + j * dfYStep, dfMinX, dfMinY,
476
0
                                dfMaxX, dfMaxY);
477
0
                        }
478
0
                        dfDXOld = dfDXNew;
479
0
                    }
480
481
0
                    dfXOld = dfX;
482
0
                    iOldOld = iOld;
483
0
                    iOld = i;
484
0
                }
485
0
            }
486
0
        }
487
0
        if (bSet)
488
0
        {
489
0
            psEnvelope->MinX = dfMinX;
490
0
            psEnvelope->MinY = dfMinY;
491
0
            psEnvelope->MaxX = dfMaxX;
492
0
            psEnvelope->MaxY = dfMaxY;
493
0
            bRet = TRUE;
494
0
        }
495
0
    }
496
497
0
    VSIFree(padfX);
498
0
    VSIFree(padfY);
499
0
    VSIFree(pabSuccess);
500
501
0
    return bRet;
502
0
}
503
504
/************************************************************************/
505
/*                             TestCapability()                         */
506
/************************************************************************/
507
508
int OGRWarpedLayer::TestCapability(const char *pszCapability)
509
0
{
510
0
    if (EQUAL(pszCapability, OLCFastGetExtent) && sStaticEnvelope.IsInit())
511
0
        return TRUE;
512
513
0
    int bVal = m_poDecoratedLayer->TestCapability(pszCapability);
514
515
0
    if (EQUAL(pszCapability, OLCFastGetArrowStream))
516
0
        return false;
517
518
0
    if (EQUAL(pszCapability, OLCFastSpatialFilter) ||
519
0
        EQUAL(pszCapability, OLCRandomWrite) ||
520
0
        EQUAL(pszCapability, OLCSequentialWrite))
521
0
    {
522
0
        if (bVal)
523
0
            bVal = m_poReversedCT != nullptr;
524
0
    }
525
0
    else if (EQUAL(pszCapability, OLCFastFeatureCount))
526
0
    {
527
0
        if (bVal)
528
0
            bVal = m_poFilterGeom == nullptr;
529
0
    }
530
531
0
    return bVal;
532
0
}
533
534
/************************************************************************/
535
/*                              SetExtent()                             */
536
/************************************************************************/
537
538
void OGRWarpedLayer::SetExtent(double dfXMin, double dfYMin, double dfXMax,
539
                               double dfYMax)
540
0
{
541
0
    sStaticEnvelope.MinX = dfXMin;
542
0
    sStaticEnvelope.MinY = dfYMin;
543
0
    sStaticEnvelope.MaxX = dfXMax;
544
0
    sStaticEnvelope.MaxY = dfYMax;
545
0
}
546
547
/************************************************************************/
548
/*                            GetArrowStream()                          */
549
/************************************************************************/
550
551
bool OGRWarpedLayer::GetArrowStream(struct ArrowArrayStream *out_stream,
552
                                    CSLConstList papszOptions)
553
0
{
554
0
    return OGRLayer::GetArrowStream(out_stream, papszOptions);
555
0
}
556
557
#endif /* #ifndef DOXYGEN_SKIP */