Coverage Report

Created: 2025-11-16 06:25

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