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 
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 
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 
78 void OGRWarpedLayer::SetSpatialFilter( OGRGeometry * poGeom )
79 {
80     SetSpatialFilter( 0, poGeom );
81 }
82 
83 /************************************************************************/
84 /*                          SetSpatialFilterRect()                      */
85 /************************************************************************/
86 
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 
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             {
GetName()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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
338 OGRErr OGRWarpedLayer::GetExtent(OGREnvelope *psExtent, int bForce)
339 {
340     return GetExtent(0, psExtent, bForce);
341 }
342 
343 /************************************************************************/
344 /*                              GetExtent()                             */
345 /************************************************************************/
346 
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 
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 
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 
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 
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 
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