1 : /******************************************************************************
2 : * $Id: ogrwarpedlayer.cpp 24633 2012-07-01 14:37:25Z 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, Even Rouault <even dot rouault at mines dash 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 24633 2012-07-01 14:37:25Z rouault $");
33 :
34 : /************************************************************************/
35 : /* OGRWarpedLayer() */
36 : /************************************************************************/
37 :
38 10 : OGRWarpedLayer::OGRWarpedLayer( OGRLayer* poDecoratedLayer,
39 : int bTakeOwnership,
40 : OGRCoordinateTransformation* poCT,
41 : OGRCoordinateTransformation* poReversedCT ) :
42 : OGRLayerDecorator(poDecoratedLayer,
43 : bTakeOwnership),
44 : m_poCT(poCT),
45 10 : m_poReversedCT(poReversedCT)
46 : {
47 10 : CPLAssert(poCT != NULL);
48 :
49 10 : if( m_poCT->GetTargetCS() != NULL )
50 : {
51 10 : m_poSRS = m_poCT->GetTargetCS();
52 10 : m_poSRS->Reference();
53 : }
54 : else
55 0 : m_poSRS = NULL;
56 10 : }
57 :
58 : /************************************************************************/
59 : /* ~OGRWarpedLayer() */
60 : /************************************************************************/
61 :
62 10 : OGRWarpedLayer::~OGRWarpedLayer()
63 : {
64 10 : if( m_poSRS != NULL )
65 10 : m_poSRS->Release();
66 10 : delete m_poCT;
67 10 : delete m_poReversedCT;
68 10 : }
69 :
70 : /************************************************************************/
71 : /* SetSpatialFilter() */
72 : /************************************************************************/
73 :
74 27 : void OGRWarpedLayer::SetSpatialFilter( OGRGeometry * poGeom )
75 : {
76 27 : OGRLayer::SetSpatialFilter(poGeom);
77 :
78 46 : if( poGeom == NULL || m_poReversedCT == NULL )
79 : {
80 19 : m_poDecoratedLayer->SetSpatialFilter(NULL);
81 : }
82 : else
83 : {
84 8 : OGREnvelope sEnvelope;
85 8 : poGeom->getEnvelope(&sEnvelope);
86 8 : if( ReprojectEnvelope(&sEnvelope, m_poReversedCT) )
87 : {
88 : m_poDecoratedLayer->SetSpatialFilterRect(sEnvelope.MinX,
89 : sEnvelope.MinY,
90 : sEnvelope.MaxX,
91 8 : sEnvelope.MaxY);
92 : }
93 : else
94 : {
95 0 : m_poDecoratedLayer->SetSpatialFilter(NULL);
96 : }
97 : }
98 27 : }
99 :
100 : /************************************************************************/
101 : /* SetSpatialFilterRect() */
102 : /************************************************************************/
103 :
104 1 : void OGRWarpedLayer::SetSpatialFilterRect( double dfMinX, double dfMinY,
105 : double dfMaxX, double dfMaxY )
106 : {
107 1 : OGRLayer::SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
108 1 : }
109 :
110 : /************************************************************************/
111 : /* GetNextFeature() */
112 : /************************************************************************/
113 :
114 300 : OGRFeature *OGRWarpedLayer::GetNextFeature()
115 : {
116 0 : while(TRUE)
117 : {
118 300 : OGRFeature* poFeature = m_poDecoratedLayer->GetNextFeature();
119 300 : if( poFeature == NULL )
120 17 : return NULL;
121 :
122 283 : OGRGeometry* poGeom = poFeature->GetGeometryRef();
123 283 : if( poGeom == NULL )
124 27 : return poFeature;
125 :
126 256 : if( poGeom->transform(m_poCT) != OGRERR_NONE )
127 : {
128 0 : delete poFeature->StealGeometry();
129 : }
130 : else
131 : {
132 256 : if( m_poFilterGeom != NULL && !FilterGeometry( poGeom ) )
133 : {
134 0 : delete poFeature;
135 0 : continue;
136 : }
137 : }
138 256 : return poFeature;
139 : }
140 : }
141 :
142 : /************************************************************************/
143 : /* GetFeature() */
144 : /************************************************************************/
145 :
146 4 : OGRFeature *OGRWarpedLayer::GetFeature( long nFID )
147 : {
148 4 : OGRFeature* poFeature = m_poDecoratedLayer->GetFeature(nFID);
149 4 : if( poFeature != NULL )
150 : {
151 2 : OGRGeometry* poGeom = poFeature->GetGeometryRef();
152 2 : if( poGeom != NULL && poGeom->transform(m_poCT) != OGRERR_NONE )
153 : {
154 0 : delete poFeature->StealGeometry();
155 : }
156 : }
157 4 : return poFeature;
158 : }
159 :
160 : /************************************************************************/
161 : /* SetFeature() */
162 : /************************************************************************/
163 :
164 1 : OGRErr OGRWarpedLayer::SetFeature( OGRFeature *poFeature )
165 : {
166 1 : OGRGeometry* poOrigGeom = poFeature->GetGeometryRef();
167 : OGRErr eErr;
168 :
169 1 : if( poOrigGeom != NULL )
170 : {
171 1 : if( m_poReversedCT == NULL )
172 0 : return OGRERR_FAILURE;
173 :
174 1 : OGRGeometry* poTransformedGeom = poOrigGeom->clone();
175 1 : if( poTransformedGeom->transform(m_poReversedCT) != OGRERR_NONE )
176 : {
177 0 : delete poTransformedGeom;
178 0 : return OGRERR_FAILURE;
179 : }
180 1 : poFeature->StealGeometry();
181 1 : poFeature->SetGeometryDirectly(poTransformedGeom);
182 :
183 1 : eErr = m_poDecoratedLayer->SetFeature(poFeature);
184 :
185 1 : poFeature->StealGeometry();
186 1 : poFeature->SetGeometryDirectly(poOrigGeom);
187 1 : delete poTransformedGeom;
188 : }
189 : else
190 0 : eErr = m_poDecoratedLayer->SetFeature(poFeature);
191 :
192 1 : return eErr;
193 : }
194 :
195 : /************************************************************************/
196 : /* CreateFeature() */
197 : /************************************************************************/
198 :
199 0 : OGRErr OGRWarpedLayer::CreateFeature( OGRFeature *poFeature )
200 : {
201 0 : OGRGeometry* poOrigGeom = poFeature->GetGeometryRef();
202 : OGRErr eErr;
203 :
204 0 : if( poOrigGeom != NULL )
205 : {
206 0 : if( m_poReversedCT == NULL )
207 0 : return OGRERR_FAILURE;
208 :
209 0 : OGRGeometry* poTransformedGeom = poOrigGeom->clone();
210 0 : if( poTransformedGeom->transform(m_poReversedCT) != OGRERR_NONE )
211 : {
212 0 : delete poTransformedGeom;
213 0 : return OGRERR_FAILURE;
214 : }
215 0 : poFeature->StealGeometry();
216 0 : poFeature->SetGeometryDirectly(poTransformedGeom);
217 :
218 0 : eErr = m_poDecoratedLayer->CreateFeature(poFeature);
219 :
220 0 : poFeature->StealGeometry();
221 0 : poFeature->SetGeometryDirectly(poOrigGeom);
222 0 : delete poTransformedGeom;
223 : }
224 : else
225 0 : eErr = m_poDecoratedLayer->CreateFeature(poFeature);
226 :
227 0 : return eErr;
228 : }
229 :
230 : /************************************************************************/
231 : /* GetSpatialRef() */
232 : /************************************************************************/
233 :
234 2 : OGRSpatialReference *OGRWarpedLayer::GetSpatialRef()
235 : {
236 2 : return m_poSRS;
237 : }
238 : /************************************************************************/
239 : /* GetFeatureCount() */
240 : /************************************************************************/
241 :
242 12 : int OGRWarpedLayer::GetFeatureCount( int bForce )
243 : {
244 12 : if( m_poFilterGeom == NULL )
245 8 : return m_poDecoratedLayer->GetFeatureCount(bForce);
246 :
247 4 : return OGRLayer::GetFeatureCount(bForce);
248 : }
249 :
250 : /************************************************************************/
251 : /* GetExtent() */
252 : /************************************************************************/
253 :
254 2 : OGRErr OGRWarpedLayer::GetExtent(OGREnvelope *psExtent, int bForce)
255 : {
256 2 : if( sStaticEnvelope.IsInit() )
257 : {
258 0 : memcpy(psExtent, &sStaticEnvelope, sizeof(OGREnvelope));
259 0 : return OGRERR_NONE;
260 : }
261 :
262 2 : OGREnvelope sExtent;
263 2 : OGRErr eErr = m_poDecoratedLayer->GetExtent(&sExtent, bForce);
264 2 : if( eErr != OGRERR_NONE )
265 0 : return eErr;
266 :
267 2 : if( ReprojectEnvelope(&sExtent, m_poCT) )
268 : {
269 2 : memcpy(psExtent, &sExtent, sizeof(OGREnvelope));
270 2 : return OGRERR_NONE;
271 : }
272 : else
273 0 : return OGRERR_FAILURE;
274 : }
275 :
276 : /************************************************************************/
277 : /* TransformAndUpdateBBAndReturnX() */
278 : /************************************************************************/
279 :
280 0 : static double TransformAndUpdateBBAndReturnX(
281 : OGRCoordinateTransformation* poCT,
282 : double dfX, double dfY,
283 : double& dfMinX, double& dfMinY, double& dfMaxX, double& dfMaxY)
284 : {
285 : int bSuccess;
286 0 : poCT->TransformEx( 1, &dfX, &dfY, NULL, &bSuccess );
287 0 : if( bSuccess )
288 : {
289 0 : if( dfX < dfMinX ) dfMinX = dfX;
290 0 : if( dfY < dfMinY ) dfMinY = dfY;
291 0 : if( dfX > dfMaxX ) dfMaxX = dfX;
292 0 : if( dfY > dfMaxY ) dfMaxY = dfY;
293 0 : return dfX;
294 : }
295 : else
296 0 : return 0.0;
297 : }
298 :
299 : /************************************************************************/
300 : /* FindXDiscontinuity() */
301 : /************************************************************************/
302 :
303 0 : static void FindXDiscontinuity(OGRCoordinateTransformation* poCT,
304 : double dfX1, double dfX2, double dfY,
305 : double& dfMinX, double& dfMinY, double& dfMaxX, double& dfMaxY,
306 : int nRecLevel = 0)
307 : {
308 0 : double dfXMid = (dfX1 + dfX2) / 2;
309 :
310 0 : double dfWrkX1 = TransformAndUpdateBBAndReturnX(poCT, dfX1, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY);
311 0 : double dfWrkXMid = TransformAndUpdateBBAndReturnX(poCT, dfXMid, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY);
312 0 : double dfWrkX2 = TransformAndUpdateBBAndReturnX(poCT, dfX2, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY);
313 :
314 0 : double dfDX1 = dfWrkXMid - dfWrkX1;
315 0 : double dfDX2 = dfWrkX2 - dfWrkXMid;
316 :
317 0 : if( dfDX1 * dfDX2 < 0 && nRecLevel < 30)
318 : {
319 0 : FindXDiscontinuity(poCT, dfX1, dfXMid, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY, nRecLevel + 1);
320 0 : FindXDiscontinuity(poCT, dfXMid, dfX2, dfY, dfMinX, dfMinY, dfMaxX, dfMaxY, nRecLevel + 1);
321 : }
322 0 : }
323 :
324 : /************************************************************************/
325 : /* ReprojectEnvelope() */
326 : /************************************************************************/
327 :
328 10 : int OGRWarpedLayer::ReprojectEnvelope( OGREnvelope* psEnvelope,
329 : OGRCoordinateTransformation* poCT )
330 : {
331 : #define NSTEP 20
332 10 : double dfXStep = (psEnvelope->MaxX - psEnvelope->MinX) / NSTEP;
333 10 : double dfYStep = (psEnvelope->MaxY - psEnvelope->MinY) / NSTEP;
334 : int i, j;
335 : double *padfX, *padfY;
336 : int* pabSuccess;
337 :
338 10 : padfX = (double*) VSIMalloc((NSTEP + 1) * (NSTEP + 1) * sizeof(double));
339 10 : padfY = (double*) VSIMalloc((NSTEP + 1) * (NSTEP + 1) * sizeof(double));
340 10 : pabSuccess = (int*) VSIMalloc((NSTEP + 1) * (NSTEP + 1) * sizeof(int));
341 10 : if( padfX == NULL || padfY == NULL || pabSuccess == NULL)
342 : {
343 0 : VSIFree(padfX);
344 0 : VSIFree(padfY);
345 0 : VSIFree(pabSuccess);
346 0 : return FALSE;
347 : }
348 :
349 220 : for(j = 0; j <= NSTEP; j++)
350 : {
351 4620 : for(i = 0; i <= NSTEP; i++)
352 : {
353 4410 : padfX[j * (NSTEP + 1) + i] = psEnvelope->MinX + i * dfXStep;
354 4410 : padfY[j * (NSTEP + 1) + i] = psEnvelope->MinY + j * dfYStep;
355 : }
356 : }
357 :
358 10 : int bRet = FALSE;
359 :
360 10 : if( poCT->TransformEx( (NSTEP + 1) * (NSTEP + 1), padfX, padfY, NULL,
361 10 : pabSuccess ) )
362 : {
363 10 : double dfMinX = 0.0, dfMinY = 0.0, dfMaxX = 0.0, dfMaxY = 0.0;
364 10 : int bSet = FALSE;
365 220 : for(j = 0; j <= NSTEP; j++)
366 : {
367 210 : double dfXOld = 0.0;
368 210 : double dfDXOld = 0.0;
369 210 : int iOld = -1, iOldOld = -1;
370 4620 : for(i = 0; i <= NSTEP; i++)
371 : {
372 4410 : if( pabSuccess[j * (NSTEP + 1) + i] )
373 : {
374 4410 : double dfX = padfX[j * (NSTEP + 1) + i];
375 4410 : double dfY = padfY[j * (NSTEP + 1) + i];
376 :
377 4410 : if( !bSet )
378 : {
379 10 : dfMinX = dfMaxX = dfX;
380 10 : dfMinY = dfMaxY = dfY;
381 10 : bSet = TRUE;
382 : }
383 : else
384 : {
385 4400 : if( dfX < dfMinX ) dfMinX = dfX;
386 4400 : if( dfY < dfMinY ) dfMinY = dfY;
387 4400 : if( dfX > dfMaxX ) dfMaxX = dfX;
388 4400 : if( dfY > dfMaxY ) dfMaxY = dfY;
389 : }
390 :
391 4410 : if( iOld >= 0 )
392 : {
393 4200 : double dfDXNew = dfX - dfXOld;
394 4200 : if( iOldOld >= 0 && dfDXNew * dfDXOld < 0 )
395 : {
396 : FindXDiscontinuity(poCT,
397 : psEnvelope->MinX + iOldOld * dfXStep,
398 : psEnvelope->MinX + i * dfXStep,
399 : psEnvelope->MinY + j * dfYStep,
400 0 : dfMinX, dfMinY, dfMaxX, dfMaxY);
401 : }
402 4200 : dfDXOld = dfDXNew;
403 : }
404 :
405 4410 : dfXOld = dfX;
406 4410 : iOldOld = iOld;
407 4410 : iOld = i;
408 :
409 : }
410 : }
411 : }
412 10 : if( bSet )
413 : {
414 10 : psEnvelope->MinX = dfMinX;
415 10 : psEnvelope->MinY = dfMinY;
416 10 : psEnvelope->MaxX = dfMaxX;
417 10 : psEnvelope->MaxY = dfMaxY;
418 10 : bRet = TRUE;
419 : }
420 : }
421 :
422 10 : VSIFree(padfX);
423 10 : VSIFree(padfY);
424 10 : VSIFree(pabSuccess);
425 :
426 10 : return bRet;
427 : }
428 :
429 : /************************************************************************/
430 : /* TestCapability() */
431 : /************************************************************************/
432 :
433 28 : int OGRWarpedLayer::TestCapability( const char * pszCapability )
434 : {
435 28 : if( EQUAL(pszCapability, OLCFastGetExtent) &&
436 : sStaticEnvelope.IsInit() )
437 0 : return TRUE;
438 :
439 28 : int bVal = m_poDecoratedLayer->TestCapability(pszCapability);
440 :
441 36 : if( EQUAL(pszCapability, OLCFastSpatialFilter) ||
442 : EQUAL(pszCapability, OLCRandomWrite) ||
443 : EQUAL(pszCapability, OLCSequentialWrite) )
444 : {
445 8 : if( bVal )
446 6 : bVal = m_poReversedCT != NULL;
447 : }
448 20 : else if( EQUAL(pszCapability, OLCFastFeatureCount) )
449 : {
450 2 : if( bVal )
451 2 : bVal = m_poFilterGeom == NULL;
452 : }
453 :
454 28 : return bVal;
455 : }
456 :
457 : /************************************************************************/
458 : /* SetExtent() */
459 : /************************************************************************/
460 :
461 0 : void OGRWarpedLayer::SetExtent( double dfXMin, double dfYMin,
462 : double dfXMax, double dfYMax )
463 : {
464 0 : sStaticEnvelope.MinX = dfXMin;
465 0 : sStaticEnvelope.MinY = dfYMin;
466 0 : sStaticEnvelope.MaxX = dfXMax;
467 0 : sStaticEnvelope.MaxY = dfYMax;
468 0 : }
|
