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