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