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