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