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 * Permission is hereby granted, free of charge, to any person obtaining a
11 * copy of this software and associated documentation files (the "Software"),
12 * to deal in the Software without restriction, including without limitation
13 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 * and/or sell copies of the Software, and to permit persons to whom the
15 * Software is furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included
18 * in all copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 * DEALINGS IN THE SOFTWARE.
27 ****************************************************************************/
28
29 #ifndef DOXYGEN_SKIP
30
31 #include "ogrwarpedlayer.h"
32
33 CPL_CVSID("$Id: ogrwarpedlayer.cpp 355b41831cd2685c85d1aabe5b95665a2c6e99b7 2019-06-19 17:07:04 +0200 Even Rouault $")
34
35 /************************************************************************/
36 /* OGRWarpedLayer() */
37 /************************************************************************/
38
OGRWarpedLayer(OGRLayer * poDecoratedLayer,int iGeomField,int bTakeOwnership,OGRCoordinateTransformation * poCT,OGRCoordinateTransformation * poReversedCT)39 OGRWarpedLayer::OGRWarpedLayer( OGRLayer* poDecoratedLayer,
40 int iGeomField,
41 int bTakeOwnership,
42 OGRCoordinateTransformation* poCT,
43 OGRCoordinateTransformation* poReversedCT ) :
44 OGRLayerDecorator(poDecoratedLayer, bTakeOwnership),
45 m_poFeatureDefn(nullptr),
46 m_iGeomField(iGeomField),
47 m_poCT(poCT),
48 m_poReversedCT(poReversedCT),
49 m_poSRS(m_poCT->GetTargetCS())
50 {
51 CPLAssert(poCT != nullptr);
52 SetDescription( poDecoratedLayer->GetDescription() );
53
54 if( m_poSRS != nullptr )
55 {
56 m_poSRS->Reference();
57 }
58 }
59
60 /************************************************************************/
61 /* ~OGRWarpedLayer() */
62 /************************************************************************/
63
~OGRWarpedLayer()64 OGRWarpedLayer::~OGRWarpedLayer()
65 {
66 if( m_poFeatureDefn != nullptr )
67 m_poFeatureDefn->Release();
68 if( m_poSRS != nullptr )
69 m_poSRS->Release();
70 delete m_poCT;
71 delete m_poReversedCT;
72 }
73
74 /************************************************************************/
75 /* SetSpatialFilter() */
76 /************************************************************************/
77
SetSpatialFilter(OGRGeometry * poGeom)78 void OGRWarpedLayer::SetSpatialFilter( OGRGeometry * poGeom )
79 {
80 SetSpatialFilter( 0, poGeom );
81 }
82
83 /************************************************************************/
84 /* SetSpatialFilterRect() */
85 /************************************************************************/
86
SetSpatialFilterRect(double dfMinX,double dfMinY,double dfMaxX,double dfMaxY)87 void OGRWarpedLayer::SetSpatialFilterRect( double dfMinX, double dfMinY,
88 double dfMaxX, double dfMaxY )
89 {
90 OGRLayer::SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
91 }
92
93 /************************************************************************/
94 /* SetSpatialFilter() */
95 /************************************************************************/
96
SetSpatialFilter(int iGeomField,OGRGeometry * poGeom)97 void OGRWarpedLayer::SetSpatialFilter( int iGeomField, OGRGeometry *poGeom )
98 {
99 if( iGeomField < 0 || iGeomField >= GetLayerDefn()->GetGeomFieldCount() )
100 {
101 CPLError(CE_Failure, CPLE_AppDefined,
102 "Invalid geometry field index : %d", iGeomField);
103 return;
104 }
105
106 m_iGeomFieldFilter = iGeomField;
107 if( InstallFilter( poGeom ) )
108 ResetReading();
109
110 if( m_iGeomFieldFilter == m_iGeomField )
111 {
112 if( poGeom == nullptr || m_poReversedCT == nullptr )
113 {
114 m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter,
115 nullptr);
116 }
117 else
118 {
119 OGREnvelope sEnvelope;
120 poGeom->getEnvelope(&sEnvelope);
121 if( CPLIsInf(sEnvelope.MinX) && CPLIsInf(sEnvelope.MinY) &&
122 CPLIsInf(sEnvelope.MaxX) && CPLIsInf(sEnvelope.MaxY) )
123 {
124 m_poDecoratedLayer->SetSpatialFilterRect(m_iGeomFieldFilter,
125 sEnvelope.MinX,
126 sEnvelope.MinY,
127 sEnvelope.MaxX,
128 sEnvelope.MaxY);
129 }
130 else if( ReprojectEnvelope(&sEnvelope, m_poReversedCT) )
131 {
132 m_poDecoratedLayer->SetSpatialFilterRect(m_iGeomFieldFilter,
133 sEnvelope.MinX,
134 sEnvelope.MinY,
135 sEnvelope.MaxX,
136 sEnvelope.MaxY);
137 }
138 else
139 {
140 m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter,
141 nullptr);
142 }
143 }
144 }
145 else
146 {
147 m_poDecoratedLayer->SetSpatialFilter(m_iGeomFieldFilter,
148 poGeom);
149 }
150 }
151
152 /************************************************************************/
153 /* SetSpatialFilterRect() */
154 /************************************************************************/
155
SetSpatialFilterRect(int iGeomField,double dfMinX,double dfMinY,double dfMaxX,double dfMaxY)156 void OGRWarpedLayer::SetSpatialFilterRect( int iGeomField, double dfMinX, double dfMinY,
157 double dfMaxX, double dfMaxY )
158 {
159 OGRLayer::SetSpatialFilterRect(iGeomField, dfMinX, dfMinY, dfMaxX, dfMaxY);
160 }
161
162 /************************************************************************/
163 /* SrcFeatureToWarpedFeature() */
164 /************************************************************************/
165
SrcFeatureToWarpedFeature(OGRFeature * poSrcFeature)166 OGRFeature *OGRWarpedLayer::SrcFeatureToWarpedFeature(OGRFeature* poSrcFeature)
167 {
168 OGRFeature* poFeature = new OGRFeature(GetLayerDefn());
169 poFeature->SetFrom(poSrcFeature);
170 poFeature->SetFID(poSrcFeature->GetFID());
171
172 OGRGeometry* poGeom = poFeature->GetGeomFieldRef(m_iGeomField);
173 if( poGeom == nullptr )
174 return poFeature;
175
176 if( poGeom->transform(m_poCT) != OGRERR_NONE )
177 {
178 delete poFeature->StealGeometry(m_iGeomField);
179 }
180
181 return poFeature;
182 }
183
184 /************************************************************************/
185 /* WarpedFeatureToSrcFeature() */
186 /************************************************************************/
187
WarpedFeatureToSrcFeature(OGRFeature * poFeature)188 OGRFeature *OGRWarpedLayer::WarpedFeatureToSrcFeature(OGRFeature* poFeature)
189 {
190 OGRFeature* poSrcFeature = new OGRFeature(m_poDecoratedLayer->GetLayerDefn());
191 poSrcFeature->SetFrom(poFeature);
192 poSrcFeature->SetFID(poFeature->GetFID());
193
194 OGRGeometry* poGeom = poSrcFeature->GetGeomFieldRef(m_iGeomField);
195 if( poGeom != nullptr )
196 {
197 if( m_poReversedCT == nullptr )
198 {
199 delete poSrcFeature;
200 return nullptr;
201 }
202
203 if( poGeom->transform(m_poReversedCT) != OGRERR_NONE )
204 {
205 delete poSrcFeature;
206 return nullptr;
207 }
208 }
209
210 return poSrcFeature;
211 }
212
213 /************************************************************************/
214 /* GetNextFeature() */
215 /************************************************************************/
216
GetNextFeature()217 OGRFeature *OGRWarpedLayer::GetNextFeature()
218 {
219 while( true )
220 {
221 OGRFeature* poFeature = m_poDecoratedLayer->GetNextFeature();
222 if( poFeature == nullptr )
223 return nullptr;
224
225 OGRFeature* poFeatureNew = SrcFeatureToWarpedFeature(poFeature);
226 delete poFeature;
227
228 OGRGeometry* poGeom = poFeatureNew->GetGeomFieldRef(m_iGeomField);
229 if( m_poFilterGeom != nullptr && !FilterGeometry( poGeom ) )
230 {
231 delete poFeatureNew;
232 continue;
233 }
234
235 return poFeatureNew;
236 }
237 }
238
239 /************************************************************************/
240 /* GetFeature() */
241 /************************************************************************/
242
GetFeature(GIntBig nFID)243 OGRFeature *OGRWarpedLayer::GetFeature( GIntBig nFID )
244 {
245 OGRFeature* poFeature = m_poDecoratedLayer->GetFeature(nFID);
246 if( poFeature != nullptr )
247 {
248 OGRFeature* poFeatureNew = SrcFeatureToWarpedFeature(poFeature);
249 delete poFeature;
250 poFeature = poFeatureNew;
251 }
252 return poFeature;
253 }
254
255 /************************************************************************/
256 /* ISetFeature() */
257 /************************************************************************/
258
ISetFeature(OGRFeature * poFeature)259 OGRErr OGRWarpedLayer::ISetFeature( OGRFeature *poFeature )
260 {
261 OGRErr eErr;
262
263 OGRFeature* poFeatureNew = WarpedFeatureToSrcFeature(poFeature);
264 if( poFeatureNew == nullptr )
265 return OGRERR_FAILURE;
266
267 eErr = m_poDecoratedLayer->SetFeature(poFeatureNew);
268
269 delete poFeatureNew;
270
271 return eErr;
272 }
273
274 /************************************************************************/
275 /* ICreateFeature() */
276 /************************************************************************/
277
ICreateFeature(OGRFeature * poFeature)278 OGRErr OGRWarpedLayer::ICreateFeature( OGRFeature *poFeature )
279 {
280 OGRErr eErr;
281
282 OGRFeature* poFeatureNew = WarpedFeatureToSrcFeature(poFeature);
283 if( poFeatureNew == nullptr )
284 return OGRERR_FAILURE;
285
286 eErr = m_poDecoratedLayer->CreateFeature(poFeatureNew);
287
288 delete poFeatureNew;
289
290 return eErr;
291 }
292
293 /************************************************************************/
294 /* GetLayerDefn() */
295 /************************************************************************/
296
GetLayerDefn()297 OGRFeatureDefn *OGRWarpedLayer::GetLayerDefn()
298 {
299 if( m_poFeatureDefn != nullptr )
300 return m_poFeatureDefn;
301
302 m_poFeatureDefn = m_poDecoratedLayer->GetLayerDefn()->Clone();
303 m_poFeatureDefn->Reference();
304 if( m_poFeatureDefn->GetGeomFieldCount() > 0 )
305 m_poFeatureDefn->GetGeomFieldDefn(m_iGeomField)->SetSpatialRef(m_poSRS);
306
307 return m_poFeatureDefn;
308 }
309
310 /************************************************************************/
311 /* GetSpatialRef() */
312 /************************************************************************/
313
GetSpatialRef()314 OGRSpatialReference *OGRWarpedLayer::GetSpatialRef()
315 {
316 if( m_iGeomField == 0 )
317 return m_poSRS;
318 else
319 return OGRLayer::GetSpatialRef();
320 }
321
322 /************************************************************************/
323 /* GetFeatureCount() */
324 /************************************************************************/
325
GetFeatureCount(int bForce)326 GIntBig OGRWarpedLayer::GetFeatureCount( int bForce )
327 {
328 if( m_poFilterGeom == nullptr )
329 return m_poDecoratedLayer->GetFeatureCount(bForce);
330
331 return OGRLayer::GetFeatureCount(bForce);
332 }
333
334 /************************************************************************/
335 /* GetExtent() */
336 /************************************************************************/
337
GetExtent(OGREnvelope * psExtent,int bForce)338 OGRErr OGRWarpedLayer::GetExtent(OGREnvelope *psExtent, int bForce)
339 {
340 return GetExtent(0, psExtent, bForce);
341 }
342
343 /************************************************************************/
344 /* GetExtent() */
345 /************************************************************************/
346
GetExtent(int iGeomField,OGREnvelope * psExtent,int bForce)347 OGRErr OGRWarpedLayer::GetExtent(int iGeomField, OGREnvelope *psExtent, int bForce)
348 {
349 if( iGeomField == m_iGeomField )
350 {
351 if( sStaticEnvelope.IsInit() )
352 {
353 *psExtent = sStaticEnvelope;
354 return OGRERR_NONE;
355 }
356
357 OGREnvelope sExtent;
358 OGRErr eErr = m_poDecoratedLayer->GetExtent(m_iGeomField, &sExtent, bForce);
359 if( eErr != OGRERR_NONE )
360 return eErr;
361
362 if( ReprojectEnvelope(&sExtent, m_poCT) )
363 {
364 *psExtent = sExtent;
365 return OGRERR_NONE;
366 }
367 else
368 return OGRERR_FAILURE;
369 }
370 else
371 return m_poDecoratedLayer->GetExtent(iGeomField, psExtent, bForce);
372 }
373
374 /************************************************************************/
375 /* TransformAndUpdateBBAndReturnX() */
376 /************************************************************************/
377
TransformAndUpdateBBAndReturnX(OGRCoordinateTransformation * poCT,double dfX,double dfY,double & dfMinX,double & dfMinY,double & dfMaxX,double & dfMaxY)378 static double TransformAndUpdateBBAndReturnX(
379 OGRCoordinateTransformation* poCT,
380 double dfX, double dfY,
381 double& dfMinX, double& dfMinY, double& dfMaxX, double& dfMaxY )
382 {
383 int bSuccess = FALSE;
384 poCT->Transform( 1, &dfX, &dfY, nullptr, nullptr, &bSuccess );
385 if( bSuccess )
386 {
387 if( dfX < dfMinX ) dfMinX = dfX;
388 if( dfY < dfMinY ) dfMinY = dfY;
389 if( dfX > dfMaxX ) dfMaxX = dfX;
390 if( dfY > dfMaxY ) dfMaxY = dfY;
391 return dfX;
392 }
393
394 return 0.0;
395 }
396
397 /************************************************************************/
398 /* FindXDiscontinuity() */
399 /************************************************************************/
400
FindXDiscontinuity(OGRCoordinateTransformation * poCT,double dfX1,double dfX2,double dfY,double & dfMinX,double & dfMinY,double & dfMaxX,double & dfMaxY,int nRecLevel=0)401 static void FindXDiscontinuity(OGRCoordinateTransformation* poCT,
402 double dfX1, double dfX2, double dfY,
403 double& dfMinX, double& dfMinY, double& dfMaxX, double& dfMaxY,
404 int nRecLevel = 0)
405 {
406 double dfXMid = (dfX1 + dfX2) / 2;
407
408 double dfWrkX1 = TransformAndUpdateBBAndReturnX(poCT, dfX1, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY);
409 double dfWrkXMid = TransformAndUpdateBBAndReturnX(poCT, dfXMid, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY);
410 double dfWrkX2 = TransformAndUpdateBBAndReturnX(poCT, dfX2, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY);
411
412 double dfDX1 = dfWrkXMid - dfWrkX1;
413 double dfDX2 = dfWrkX2 - dfWrkXMid;
414
415 if( dfDX1 * dfDX2 < 0 && nRecLevel < 30)
416 {
417 FindXDiscontinuity(poCT, dfX1, dfXMid, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY, nRecLevel + 1);
418 FindXDiscontinuity(poCT, dfXMid, dfX2, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY, nRecLevel + 1);
419 }
420 }
421
422 /************************************************************************/
423 /* ReprojectEnvelope() */
424 /************************************************************************/
425
ReprojectEnvelope(OGREnvelope * psEnvelope,OGRCoordinateTransformation * poCT)426 int OGRWarpedLayer::ReprojectEnvelope( OGREnvelope* psEnvelope,
427 OGRCoordinateTransformation* poCT )
428 {
429 const int NSTEP = 20;
430 double dfXStep = (psEnvelope->MaxX - psEnvelope->MinX) / NSTEP;
431 double dfYStep = (psEnvelope->MaxY - psEnvelope->MinY) / NSTEP;
432
433 double *padfX = static_cast<double*>(
434 VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(double)));
435 double *padfY = static_cast<double*>(
436 VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(double)));
437 int* pabSuccess = static_cast<int*>(
438 VSI_MALLOC_VERBOSE((NSTEP + 1) * (NSTEP + 1) * sizeof(int)));
439 if( padfX == nullptr || padfY == nullptr || pabSuccess == nullptr)
440 {
441 VSIFree(padfX);
442 VSIFree(padfY);
443 VSIFree(pabSuccess);
444 return FALSE;
445 }
446
447 for( int j = 0; j <= NSTEP; j++ )
448 {
449 for( int i = 0; i <= NSTEP; i++ )
450 {
451 padfX[j * (NSTEP + 1) + i] = psEnvelope->MinX + i * dfXStep;
452 padfY[j * (NSTEP + 1) + i] = psEnvelope->MinY + j * dfYStep;
453 }
454 }
455
456 int bRet = FALSE;
457
458 if( poCT->Transform( (NSTEP + 1) * (NSTEP + 1), padfX, padfY, nullptr,
459 nullptr, pabSuccess ) )
460 {
461 double dfMinX = 0.0;
462 double dfMinY = 0.0;
463 double dfMaxX = 0.0;
464 double dfMaxY = 0.0;
465 int bSet = FALSE;
466 for( int j = 0; j <= NSTEP; j++ )
467 {
468 double dfXOld = 0.0;
469 double dfDXOld = 0.0;
470 int iOld = -1;
471 int iOldOld = -1;
472 for( int i = 0; i <= NSTEP; i++ )
473 {
474 if( pabSuccess[j * (NSTEP + 1) + i] )
475 {
476 double dfX = padfX[j * (NSTEP + 1) + i];
477 double dfY = padfY[j * (NSTEP + 1) + i];
478
479 if( !bSet )
480 {
481 dfMinX = dfX;
482 dfMaxX = dfX;
483 dfMinY = dfY;
484 dfMaxY = dfY;
485 bSet = TRUE;
486 }
487 else
488 {
489 if( dfX < dfMinX ) dfMinX = dfX;
490 if( dfY < dfMinY ) dfMinY = dfY;
491 if( dfX > dfMaxX ) dfMaxX = dfX;
492 if( dfY > dfMaxY ) dfMaxY = dfY;
493 }
494
495 if( iOld >= 0 )
496 {
497 double dfDXNew = dfX - dfXOld;
498 if( iOldOld >= 0 && dfDXNew * dfDXOld < 0 )
499 {
500 FindXDiscontinuity(poCT,
501 psEnvelope->MinX + iOldOld * dfXStep,
502 psEnvelope->MinX + i * dfXStep,
503 psEnvelope->MinY + j * dfYStep,
504 dfMinX, dfMinY, dfMaxX, dfMaxY);
505 }
506 dfDXOld = dfDXNew;
507 }
508
509 dfXOld = dfX;
510 iOldOld = iOld;
511 iOld = i;
512 }
513 }
514 }
515 if( bSet )
516 {
517 psEnvelope->MinX = dfMinX;
518 psEnvelope->MinY = dfMinY;
519 psEnvelope->MaxX = dfMaxX;
520 psEnvelope->MaxY = dfMaxY;
521 bRet = TRUE;
522 }
523 }
524
525 VSIFree(padfX);
526 VSIFree(padfY);
527 VSIFree(pabSuccess);
528
529 return bRet;
530 }
531
532 /************************************************************************/
533 /* TestCapability() */
534 /************************************************************************/
535
TestCapability(const char * pszCapability)536 int OGRWarpedLayer::TestCapability( const char * pszCapability )
537 {
538 if( EQUAL(pszCapability, OLCFastGetExtent) &&
539 sStaticEnvelope.IsInit() )
540 return TRUE;
541
542 int bVal = m_poDecoratedLayer->TestCapability(pszCapability);
543
544 if( EQUAL(pszCapability, OLCFastSpatialFilter) ||
545 EQUAL(pszCapability, OLCRandomWrite) ||
546 EQUAL(pszCapability, OLCSequentialWrite) )
547 {
548 if( bVal )
549 bVal = m_poReversedCT != nullptr;
550 }
551 else if( EQUAL(pszCapability, OLCFastFeatureCount) )
552 {
553 if( bVal )
554 bVal = m_poFilterGeom == nullptr;
555 }
556
557 return bVal;
558 }
559
560 /************************************************************************/
561 /* SetExtent() */
562 /************************************************************************/
563
SetExtent(double dfXMin,double dfYMin,double dfXMax,double dfYMax)564 void OGRWarpedLayer::SetExtent( double dfXMin, double dfYMin,
565 double dfXMax, double dfYMax )
566 {
567 sStaticEnvelope.MinX = dfXMin;
568 sStaticEnvelope.MinY = dfYMin;
569 sStaticEnvelope.MaxX = dfXMax;
570 sStaticEnvelope.MaxY = dfYMax;
571 }
572
573 #endif /* #ifndef DOXYGEN_SKIP */
574