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