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