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 17 : OGRWarpedLayer::OGRWarpedLayer(OGRLayer *poDecoratedLayer, int iGeomField,
24 : int bTakeOwnership,
25 : OGRCoordinateTransformation *poCT,
26 17 : OGRCoordinateTransformation *poReversedCT)
27 : : OGRLayerDecorator(poDecoratedLayer, bTakeOwnership),
28 : m_poFeatureDefn(nullptr), m_iGeomField(iGeomField), m_poCT(poCT),
29 : m_poReversedCT(poReversedCT),
30 17 : m_poSRS(const_cast<OGRSpatialReference *>(m_poCT->GetTargetCS()))
31 : {
32 17 : CPLAssert(poCT != nullptr);
33 17 : SetDescription(poDecoratedLayer->GetDescription());
34 :
35 17 : if (m_poSRS != nullptr)
36 : {
37 17 : m_poSRS->Reference();
38 : }
39 17 : }
40 :
41 : /************************************************************************/
42 : /* ~OGRWarpedLayer() */
43 : /************************************************************************/
44 :
45 34 : OGRWarpedLayer::~OGRWarpedLayer()
46 : {
47 17 : if (m_poFeatureDefn != nullptr)
48 14 : m_poFeatureDefn->Release();
49 17 : if (m_poSRS != nullptr)
50 17 : m_poSRS->Release();
51 17 : delete m_poCT;
52 17 : delete m_poReversedCT;
53 34 : }
54 :
55 : /************************************************************************/
56 : /* SetSpatialFilter() */
57 : /************************************************************************/
58 :
59 15 : void OGRWarpedLayer::SetSpatialFilter(OGRGeometry *poGeom)
60 : {
61 15 : SetSpatialFilter(0, poGeom);
62 15 : }
63 :
64 : /************************************************************************/
65 : /* SetSpatialFilterRect() */
66 : /************************************************************************/
67 :
68 1 : void OGRWarpedLayer::SetSpatialFilterRect(double dfMinX, double dfMinY,
69 : double dfMaxX, double dfMaxY)
70 : {
71 1 : OGRLayer::SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
72 1 : }
73 :
74 : /************************************************************************/
75 : /* SetSpatialFilter() */
76 : /************************************************************************/
77 :
78 221 : void OGRWarpedLayer::SetSpatialFilter(int iGeomField, OGRGeometry *poGeom)
79 : {
80 221 : if (iGeomField < 0 || iGeomField >= GetLayerDefn()->GetGeomFieldCount())
81 : {
82 2 : CPLError(CE_Failure, CPLE_AppDefined,
83 : "Invalid geometry field index : %d", iGeomField);
84 2 : return;
85 : }
86 :
87 219 : m_iGeomFieldFilter = iGeomField;
88 219 : if (InstallFilter(poGeom))
89 99 : ResetReading();
90 :
91 219 : if (m_iGeomFieldFilter == m_iGeomField)
92 : {
93 63 : if (poGeom == nullptr || m_poReversedCT == nullptr)
94 : {
95 19 : m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter, nullptr);
96 : }
97 : else
98 : {
99 44 : OGREnvelope sEnvelope;
100 44 : poGeom->getEnvelope(&sEnvelope);
101 54 : if (std::isinf(sEnvelope.MinX) && std::isinf(sEnvelope.MinY) &&
102 54 : std::isinf(sEnvelope.MaxX) && std::isinf(sEnvelope.MaxY))
103 : {
104 5 : m_poDecoratedLayer->SetSpatialFilterRect(
105 : m_iGeomFieldFilter, sEnvelope.MinX, sEnvelope.MinY,
106 5 : sEnvelope.MaxX, sEnvelope.MaxY);
107 : }
108 39 : else if (ReprojectEnvelope(&sEnvelope, m_poReversedCT))
109 : {
110 39 : m_poDecoratedLayer->SetSpatialFilterRect(
111 : m_iGeomFieldFilter, sEnvelope.MinX, sEnvelope.MinY,
112 39 : sEnvelope.MaxX, sEnvelope.MaxY);
113 : }
114 : else
115 : {
116 0 : m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter,
117 0 : nullptr);
118 : }
119 : }
120 : }
121 : else
122 : {
123 156 : m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter, poGeom);
124 : }
125 : }
126 :
127 : /************************************************************************/
128 : /* SetSpatialFilterRect() */
129 : /************************************************************************/
130 :
131 1 : void OGRWarpedLayer::SetSpatialFilterRect(int iGeomField, double dfMinX,
132 : double dfMinY, double dfMaxX,
133 : double dfMaxY)
134 : {
135 1 : OGRLayer::SetSpatialFilterRect(iGeomField, dfMinX, dfMinY, dfMaxX, dfMaxY);
136 1 : }
137 :
138 : /************************************************************************/
139 : /* SrcFeatureToWarpedFeature() */
140 : /************************************************************************/
141 :
142 : std::unique_ptr<OGRFeature>
143 398 : OGRWarpedLayer::SrcFeatureToWarpedFeature(std::unique_ptr<OGRFeature> poFeature)
144 : {
145 : // This is safe to do here as they have matching attribute and geometry
146 : // fields
147 398 : poFeature->SetFDefnUnsafe(GetLayerDefn());
148 :
149 398 : OGRGeometry *poGeom = poFeature->GetGeomFieldRef(m_iGeomField);
150 398 : if (poGeom && poGeom->transform(m_poCT) != OGRERR_NONE)
151 : {
152 2 : delete poFeature->StealGeometry(m_iGeomField);
153 : }
154 :
155 398 : return poFeature;
156 : }
157 :
158 : /************************************************************************/
159 : /* WarpedFeatureToSrcFeature() */
160 : /************************************************************************/
161 :
162 : std::unique_ptr<OGRFeature>
163 9 : OGRWarpedLayer::WarpedFeatureToSrcFeature(std::unique_ptr<OGRFeature> poFeature)
164 : {
165 : // This is safe to do here as they have matching attribute and geometry
166 : // fields
167 9 : poFeature->SetFDefnUnsafe(m_poDecoratedLayer->GetLayerDefn());
168 :
169 9 : OGRGeometry *poGeom = poFeature->GetGeomFieldRef(m_iGeomField);
170 16 : if (poGeom &&
171 7 : (!m_poReversedCT || poGeom->transform(m_poReversedCT) != OGRERR_NONE))
172 : {
173 2 : return nullptr;
174 : }
175 :
176 7 : return poFeature;
177 : }
178 :
179 : /************************************************************************/
180 : /* GetNextFeature() */
181 : /************************************************************************/
182 :
183 487 : OGRFeature *OGRWarpedLayer::GetNextFeature()
184 : {
185 : while (true)
186 : {
187 : auto poFeature =
188 487 : std::unique_ptr<OGRFeature>(m_poDecoratedLayer->GetNextFeature());
189 487 : if (!poFeature)
190 95 : return nullptr;
191 :
192 392 : auto poFeatureNew = SrcFeatureToWarpedFeature(std::move(poFeature));
193 392 : const OGRGeometry *poGeom = poFeatureNew->GetGeomFieldRef(m_iGeomField);
194 392 : if (m_poFilterGeom != nullptr && !FilterGeometry(poGeom))
195 : {
196 0 : continue;
197 : }
198 :
199 392 : return poFeatureNew.release();
200 0 : }
201 : }
202 :
203 : /************************************************************************/
204 : /* GetFeature() */
205 : /************************************************************************/
206 :
207 9 : OGRFeature *OGRWarpedLayer::GetFeature(GIntBig nFID)
208 : {
209 : auto poFeature =
210 18 : std::unique_ptr<OGRFeature>(m_poDecoratedLayer->GetFeature(nFID));
211 9 : if (poFeature)
212 : {
213 6 : poFeature = SrcFeatureToWarpedFeature(std::move(poFeature));
214 : }
215 18 : return poFeature.release();
216 : }
217 :
218 : /************************************************************************/
219 : /* ISetFeature() */
220 : /************************************************************************/
221 :
222 5 : OGRErr OGRWarpedLayer::ISetFeature(OGRFeature *poFeature)
223 : {
224 : auto poFeatureNew = WarpedFeatureToSrcFeature(
225 10 : std::unique_ptr<OGRFeature>(poFeature->Clone()));
226 5 : if (!poFeatureNew)
227 1 : return OGRERR_FAILURE;
228 :
229 4 : return m_poDecoratedLayer->SetFeature(poFeatureNew.get());
230 : }
231 :
232 : /************************************************************************/
233 : /* ICreateFeature() */
234 : /************************************************************************/
235 :
236 4 : OGRErr OGRWarpedLayer::ICreateFeature(OGRFeature *poFeature)
237 : {
238 : auto poFeatureNew = WarpedFeatureToSrcFeature(
239 8 : std::unique_ptr<OGRFeature>(poFeature->Clone()));
240 4 : if (!poFeatureNew)
241 1 : return OGRERR_FAILURE;
242 :
243 3 : return m_poDecoratedLayer->CreateFeature(poFeatureNew.get());
244 : }
245 :
246 : /************************************************************************/
247 : /* IUpsertFeature() */
248 : /************************************************************************/
249 :
250 0 : OGRErr OGRWarpedLayer::IUpsertFeature(OGRFeature *poFeature)
251 : {
252 : auto poFeatureNew = WarpedFeatureToSrcFeature(
253 0 : std::unique_ptr<OGRFeature>(poFeature->Clone()));
254 0 : if (poFeatureNew == nullptr)
255 0 : return OGRERR_FAILURE;
256 :
257 0 : return m_poDecoratedLayer->UpsertFeature(poFeatureNew.get());
258 : }
259 :
260 : /************************************************************************/
261 : /* IUpdateFeature() */
262 : /************************************************************************/
263 :
264 0 : OGRErr OGRWarpedLayer::IUpdateFeature(OGRFeature *poFeature,
265 : int nUpdatedFieldsCount,
266 : const int *panUpdatedFieldsIdx,
267 : int nUpdatedGeomFieldsCount,
268 : const int *panUpdatedGeomFieldsIdx,
269 : bool bUpdateStyleString)
270 : {
271 : auto poFeatureNew = WarpedFeatureToSrcFeature(
272 0 : std::unique_ptr<OGRFeature>(poFeature->Clone()));
273 0 : if (!poFeatureNew)
274 0 : return OGRERR_FAILURE;
275 :
276 0 : return m_poDecoratedLayer->UpdateFeature(
277 : poFeatureNew.get(), nUpdatedFieldsCount, panUpdatedFieldsIdx,
278 0 : nUpdatedGeomFieldsCount, panUpdatedGeomFieldsIdx, bUpdateStyleString);
279 : }
280 :
281 : /************************************************************************/
282 : /* GetLayerDefn() */
283 : /************************************************************************/
284 :
285 1100 : OGRFeatureDefn *OGRWarpedLayer::GetLayerDefn()
286 : {
287 1100 : if (m_poFeatureDefn != nullptr)
288 1086 : return m_poFeatureDefn;
289 :
290 14 : m_poFeatureDefn = m_poDecoratedLayer->GetLayerDefn()->Clone();
291 14 : m_poFeatureDefn->Reference();
292 14 : if (m_poFeatureDefn->GetGeomFieldCount() > 0)
293 14 : m_poFeatureDefn->GetGeomFieldDefn(m_iGeomField)->SetSpatialRef(m_poSRS);
294 :
295 14 : return m_poFeatureDefn;
296 : }
297 :
298 : /************************************************************************/
299 : /* GetSpatialRef() */
300 : /************************************************************************/
301 :
302 5 : OGRSpatialReference *OGRWarpedLayer::GetSpatialRef()
303 : {
304 5 : if (m_iGeomField == 0)
305 5 : return m_poSRS;
306 : else
307 0 : return OGRLayer::GetSpatialRef();
308 : }
309 :
310 : /************************************************************************/
311 : /* GetFeatureCount() */
312 : /************************************************************************/
313 :
314 45 : GIntBig OGRWarpedLayer::GetFeatureCount(int bForce)
315 : {
316 45 : if (m_poFilterGeom == nullptr)
317 25 : return m_poDecoratedLayer->GetFeatureCount(bForce);
318 :
319 20 : return OGRLayer::GetFeatureCount(bForce);
320 : }
321 :
322 : /************************************************************************/
323 : /* GetExtent() */
324 : /************************************************************************/
325 :
326 3 : OGRErr OGRWarpedLayer::GetExtent(OGREnvelope *psExtent, int bForce)
327 : {
328 3 : return GetExtent(0, psExtent, bForce);
329 : }
330 :
331 : /************************************************************************/
332 : /* GetExtent() */
333 : /************************************************************************/
334 :
335 13 : OGRErr OGRWarpedLayer::GetExtent(int iGeomField, OGREnvelope *psExtent,
336 : int bForce)
337 : {
338 13 : if (iGeomField == m_iGeomField)
339 : {
340 7 : if (sStaticEnvelope.IsInit())
341 : {
342 1 : *psExtent = sStaticEnvelope;
343 1 : return OGRERR_NONE;
344 : }
345 :
346 6 : OGREnvelope sExtent;
347 : OGRErr eErr =
348 6 : m_poDecoratedLayer->GetExtent(m_iGeomField, &sExtent, bForce);
349 6 : if (eErr != OGRERR_NONE)
350 0 : return eErr;
351 :
352 6 : if (ReprojectEnvelope(&sExtent, m_poCT))
353 : {
354 6 : *psExtent = sExtent;
355 6 : return OGRERR_NONE;
356 : }
357 : else
358 0 : return OGRERR_FAILURE;
359 : }
360 : else
361 6 : 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 45 : int OGRWarpedLayer::ReprojectEnvelope(OGREnvelope *psEnvelope,
426 : OGRCoordinateTransformation *poCT)
427 : {
428 45 : const int NSTEP = 20;
429 45 : double dfXStep = (psEnvelope->MaxX - psEnvelope->MinX) / NSTEP;
430 45 : double dfYStep = (psEnvelope->MaxY - psEnvelope->MinY) / NSTEP;
431 :
432 : double *padfX = static_cast<double *>(
433 45 : VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(double)));
434 : double *padfY = static_cast<double *>(
435 45 : VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(double)));
436 : int *pabSuccess = static_cast<int *>(
437 45 : VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(int)));
438 45 : 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 990 : for (int j = 0; j <= NSTEP; j++)
447 : {
448 20790 : for (int i = 0; i <= NSTEP; i++)
449 : {
450 19845 : padfX[j * (NSTEP + 1) + i] = psEnvelope->MinX + i * dfXStep;
451 19845 : padfY[j * (NSTEP + 1) + i] = psEnvelope->MinY + j * dfYStep;
452 : }
453 : }
454 :
455 45 : int bRet = FALSE;
456 :
457 90 : if (poCT->Transform((NSTEP + 1) * (NSTEP + 1), padfX, padfY, nullptr,
458 45 : nullptr, pabSuccess))
459 : {
460 45 : double dfMinX = 0.0;
461 45 : double dfMinY = 0.0;
462 45 : double dfMaxX = 0.0;
463 45 : double dfMaxY = 0.0;
464 45 : int bSet = FALSE;
465 990 : for (int j = 0; j <= NSTEP; j++)
466 : {
467 945 : double dfXOld = 0.0;
468 945 : double dfDXOld = 0.0;
469 945 : int iOld = -1;
470 945 : int iOldOld = -1;
471 20790 : for (int i = 0; i <= NSTEP; i++)
472 : {
473 19845 : if (pabSuccess[j * (NSTEP + 1) + i])
474 : {
475 17665 : double dfX = padfX[j * (NSTEP + 1) + i];
476 17665 : double dfY = padfY[j * (NSTEP + 1) + i];
477 :
478 17665 : if (!bSet)
479 : {
480 45 : dfMinX = dfX;
481 45 : dfMaxX = dfX;
482 45 : dfMinY = dfY;
483 45 : dfMaxY = dfY;
484 45 : bSet = TRUE;
485 : }
486 : else
487 : {
488 17620 : if (dfX < dfMinX)
489 124 : dfMinX = dfX;
490 17620 : if (dfY < dfMinY)
491 81 : dfMinY = dfY;
492 17620 : if (dfX > dfMaxX)
493 890 : dfMaxX = dfX;
494 17620 : if (dfY > dfMaxY)
495 3393 : dfMaxY = dfY;
496 : }
497 :
498 17665 : if (iOld >= 0)
499 : {
500 16800 : double dfDXNew = dfX - dfXOld;
501 16800 : 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 16800 : dfDXOld = dfDXNew;
510 : }
511 :
512 17665 : dfXOld = dfX;
513 17665 : iOldOld = iOld;
514 17665 : iOld = i;
515 : }
516 : }
517 : }
518 45 : if (bSet)
519 : {
520 45 : psEnvelope->MinX = dfMinX;
521 45 : psEnvelope->MinY = dfMinY;
522 45 : psEnvelope->MaxX = dfMaxX;
523 45 : psEnvelope->MaxY = dfMaxY;
524 45 : bRet = TRUE;
525 : }
526 : }
527 :
528 45 : VSIFree(padfX);
529 45 : VSIFree(padfY);
530 45 : VSIFree(pabSuccess);
531 :
532 45 : return bRet;
533 : }
534 :
535 : /************************************************************************/
536 : /* TestCapability() */
537 : /************************************************************************/
538 :
539 174 : int OGRWarpedLayer::TestCapability(const char *pszCapability)
540 : {
541 174 : if (EQUAL(pszCapability, OLCFastGetExtent) && sStaticEnvelope.IsInit())
542 0 : return TRUE;
543 :
544 174 : int bVal = m_poDecoratedLayer->TestCapability(pszCapability);
545 :
546 174 : if (EQUAL(pszCapability, OLCFastGetArrowStream))
547 1 : return false;
548 :
549 173 : if (EQUAL(pszCapability, OLCFastSpatialFilter) ||
550 173 : EQUAL(pszCapability, OLCRandomWrite) ||
551 172 : EQUAL(pszCapability, OLCSequentialWrite))
552 : {
553 2 : if (bVal)
554 0 : bVal = m_poReversedCT != nullptr;
555 : }
556 171 : else if (EQUAL(pszCapability, OLCFastFeatureCount))
557 : {
558 0 : if (bVal)
559 0 : bVal = m_poFilterGeom == nullptr;
560 : }
561 :
562 173 : 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 */
|
