Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  OpenGIS Simple Features Reference Implementation
       4             :  * Purpose:  Implementation of PMTiles
       5             :  * Author:   Even Rouault <even.rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2023, Planet Labs
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "ogr_pmtiles.h"
      14             : 
      15             : #include <algorithm>
      16             : #include <limits>
      17             : 
      18             : /************************************************************************/
      19             : /*                   find_tile_idx_lesser_or_equal()                    */
      20             : /************************************************************************/
      21             : 
      22             : static int
      23        1270 : find_tile_idx_lesser_or_equal(const std::vector<pmtiles::entryv3> &entries,
      24             :                               uint64_t tile_id)
      25             : {
      26        1270 :     if (!entries.empty() && tile_id <= entries[0].tile_id)
      27         253 :         return 0;
      28             : 
      29        1017 :     int m = 0;
      30        1017 :     int n = static_cast<int>(entries.size()) - 1;
      31        3891 :     while (m <= n)
      32             :     {
      33        2965 :         int k = (n + m) >> 1;
      34        2965 :         if (tile_id > entries[k].tile_id)
      35             :         {
      36        1912 :             m = k + 1;
      37             :         }
      38        1053 :         else if (tile_id < entries[k].tile_id)
      39             :         {
      40         962 :             n = k - 1;
      41             :         }
      42             :         else
      43             :         {
      44          91 :             return k;
      45             :         }
      46             :     }
      47             : 
      48         926 :     return n;
      49             : }
      50             : 
      51             : /************************************************************************/
      52             : /*                         LoadRootDirectory()                          */
      53             : /************************************************************************/
      54             : 
      55         379 : bool OGRPMTilesTileIterator::LoadRootDirectory()
      56             : {
      57         379 :     if (m_nZoomLevel >= 0)
      58             :     {
      59         379 :         CPLDebugOnly("PMTiles", "minx=%d miny=%d maxx=%d maxy=%d", m_nMinX,
      60             :                      m_nMinY, m_nMaxX, m_nMaxY);
      61             :         // If we don't query too many tiles, establish the minimum
      62             :         // and maximum tile id, we are interested in.
      63             :         // (is there a clever way of figuring out this?)
      64         379 :         if (m_nMinX >= 0 && m_nMinY >= 0 && m_nMaxX >= m_nMinX &&
      65         296 :             m_nMaxY >= m_nMinY &&
      66         296 :             (m_nMaxX - m_nMinX + 1) <= 100 / (m_nMaxY - m_nMinY + 1))
      67             :         {
      68         376 :             for (int iY = m_nMinY; iY <= m_nMaxY; ++iY)
      69             :             {
      70         843 :                 for (int iX = m_nMinX; iX <= m_nMaxX; ++iX)
      71             :                 {
      72        1252 :                     const uint64_t nTileId = pmtiles::zxy_to_tileid(
      73         627 :                         static_cast<uint8_t>(m_nZoomLevel), iX, iY);
      74         625 :                     m_nMinTileId = std::min(m_nMinTileId, nTileId);
      75         625 :                     m_nMaxTileId = std::max(m_nMaxTileId, nTileId);
      76             :                 }
      77         158 :             }
      78             :         }
      79             :         else
      80             :         {
      81         438 :             m_nMinTileId = pmtiles::zxy_to_tileid(
      82         219 :                 static_cast<uint8_t>(m_nZoomLevel), 0, 0);
      83         438 :             m_nMaxTileId = pmtiles::zxy_to_tileid(
      84         219 :                                static_cast<uint8_t>(m_nZoomLevel) + 1, 0, 0) -
      85             :                            1;
      86             :         }
      87             : 
      88             :         // If filtering by bbox and that the gap between min and max
      89             :         // tile id is too big, use an iteration over (x, y) space rather
      90             :         // than tile id space.
      91             : 
      92             :         // Config option for debugging purposes
      93         377 :         const unsigned nThreshold = static_cast<unsigned>(atoi(
      94         377 :             CPLGetConfigOption("OGR_PMTILES_ITERATOR_THRESHOLD", "10000")));
      95         377 :         if (m_nMinX >= 0 && m_nMinY >= 0 && m_nMaxX >= m_nMinX &&
      96         294 :             m_nMaxY >= m_nMinY &&
      97         294 :             !(m_nMinX == 0 && m_nMinY == 0 &&
      98         157 :               m_nMaxX == (1 << m_nZoomLevel) - 1 &&
      99         144 :               m_nMaxY == (1 << m_nZoomLevel) - 1) &&
     100         152 :             m_nMaxTileId - m_nMinTileId > nThreshold)
     101             :         {
     102           4 :             m_nCurX = m_nMinX;
     103           4 :             m_nCurY = m_nMinY;
     104           8 :             m_nMinTileId = pmtiles::zxy_to_tileid(
     105           4 :                 static_cast<uint8_t>(m_nZoomLevel), m_nCurX, m_nCurY);
     106           4 :             m_nMaxTileId = m_nMinTileId;
     107             :         }
     108             :     }
     109             : 
     110         377 :     const auto &sHeader = m_poDS->GetHeader();
     111         754 :     const auto *posStr = m_poDS->ReadInternal(
     112         377 :         sHeader.root_dir_offset, static_cast<uint32_t>(sHeader.root_dir_bytes),
     113             :         "header");
     114         377 :     if (!posStr)
     115             :     {
     116           0 :         return false;
     117             :     }
     118             : 
     119         754 :     DirectoryContext sContext;
     120         377 :     sContext.sEntries = pmtiles::deserialize_directory(*posStr);
     121             : 
     122         377 :     if (m_nZoomLevel >= 0)
     123             :     {
     124         377 :         if (m_nCurX >= 0)
     125             :         {
     126             :             while (true)
     127             :             {
     128           4 :                 const int nMinEntryIdx = find_tile_idx_lesser_or_equal(
     129             :                     sContext.sEntries, m_nMinTileId);
     130           4 :                 if (nMinEntryIdx < 0)
     131             :                 {
     132           0 :                     m_nCurX++;
     133           0 :                     if (m_nCurX > m_nMaxX)
     134             :                     {
     135           0 :                         m_nCurX = m_nMinX;
     136           0 :                         m_nCurY++;
     137           0 :                         if (m_nCurY > m_nMaxY)
     138             :                         {
     139           0 :                             return false;
     140             :                         }
     141             :                     }
     142           0 :                     m_nMinTileId = pmtiles::zxy_to_tileid(
     143           0 :                         static_cast<uint8_t>(m_nZoomLevel), m_nCurX, m_nCurY);
     144           0 :                     m_nMaxTileId = m_nMinTileId;
     145             :                 }
     146             :                 else
     147             :                 {
     148           4 :                     sContext.nIdxInEntries = nMinEntryIdx;
     149           4 :                     break;
     150             :                 }
     151           0 :             }
     152             :         }
     153             :         else
     154             :         {
     155             :             const int nMinEntryIdx =
     156         373 :                 find_tile_idx_lesser_or_equal(sContext.sEntries, m_nMinTileId);
     157         373 :             if (nMinEntryIdx < 0)
     158             :             {
     159           0 :                 return false;
     160             :             }
     161         373 :             sContext.nIdxInEntries = nMinEntryIdx;
     162             :         }
     163             :     }
     164             : 
     165         377 :     m_aoStack.emplace(std::move(sContext));
     166         377 :     return true;
     167             : }
     168             : 
     169             : /************************************************************************/
     170             : /*                           SkipRunLength()                            */
     171             : /************************************************************************/
     172             : 
     173           2 : void OGRPMTilesTileIterator::SkipRunLength()
     174             : {
     175           2 :     if (!m_aoStack.empty())
     176             :     {
     177           2 :         auto &topContext = m_aoStack.top();
     178           2 :         if (topContext.nIdxInEntries < topContext.sEntries.size())
     179             :         {
     180             :             const auto &sCurrentEntry =
     181           2 :                 topContext.sEntries[topContext.nIdxInEntries];
     182           2 :             if (sCurrentEntry.run_length > 1)
     183             :             {
     184           2 :                 m_nLastTileId =
     185           2 :                     sCurrentEntry.tile_id + sCurrentEntry.run_length - 1;
     186           2 :                 topContext.nIdxInRunLength = sCurrentEntry.run_length;
     187             :             }
     188             :         }
     189             :     }
     190           2 : }
     191             : 
     192             : /************************************************************************/
     193             : /*                            GetNextTile()                             */
     194             : /************************************************************************/
     195             : 
     196        1057 : pmtiles::entry_zxy OGRPMTilesTileIterator::GetNextTile(uint32_t *pnRunLength)
     197             : {
     198        1057 :     if (m_bEOF)
     199          15 :         return pmtiles::entry_zxy(0, 0, 0, 0, 0);
     200             : 
     201        1042 :     const auto &sHeader = m_poDS->GetHeader();
     202             :     try
     203             :     {
     204             :         // Put the root directory as the first element of the stack
     205             :         // of directories, if the stack is empty
     206        1042 :         if (m_aoStack.empty())
     207             :         {
     208         379 :             if (!LoadRootDirectory())
     209             :             {
     210           0 :                 m_bEOF = true;
     211         850 :                 return pmtiles::entry_zxy(0, 0, 0, 0, 0);
     212             :             }
     213             :         }
     214             : 
     215        7171 :         const auto AdvanceToNextTile = [this]()
     216             :         {
     217        1704 :             if (m_nCurX >= 0)
     218             :             {
     219             :                 while (true)
     220             :                 {
     221         682 :                     m_nCurX++;
     222         682 :                     if (m_nCurX > m_nMaxX)
     223             :                     {
     224          35 :                         m_nCurX = m_nMinX;
     225          35 :                         m_nCurY++;
     226          35 :                         if (m_nCurY > m_nMaxY)
     227             :                         {
     228           4 :                             m_bEOF = true;
     229           4 :                             return false;
     230             :                         }
     231             :                     }
     232         678 :                     if (!m_bEOF)
     233             :                     {
     234        1356 :                         m_nMinTileId = pmtiles::zxy_to_tileid(
     235         678 :                             static_cast<uint8_t>(m_nZoomLevel), m_nCurX,
     236         678 :                             m_nCurY);
     237         678 :                         m_nMaxTileId = m_nMinTileId;
     238         678 :                         m_nLastTileId = INVALID_LAST_TILE_ID;
     239         678 :                         while (m_aoStack.size() > 1)
     240           0 :                             m_aoStack.pop();
     241        1356 :                         const int nMinEntryIdx = find_tile_idx_lesser_or_equal(
     242         678 :                             m_aoStack.top().sEntries, m_nMinTileId);
     243         678 :                         if (nMinEntryIdx < 0)
     244             :                         {
     245           0 :                             continue;
     246             :                         }
     247         678 :                         m_aoStack.top().nIdxInEntries = nMinEntryIdx;
     248         678 :                         m_aoStack.top().nIdxInRunLength = 0;
     249         678 :                         break;
     250             :                     }
     251           0 :                 }
     252         678 :                 return true;
     253             :             }
     254        1022 :             return false;
     255        1040 :         };
     256             : 
     257             :         while (true)
     258             :         {
     259        9120 :             if (m_aoStack.top().nIdxInEntries ==
     260        4560 :                 m_aoStack.top().sEntries.size())
     261             :             {
     262         433 :                 if (m_aoStack.size() == 1 && AdvanceToNextTile())
     263         129 :                     continue;
     264             : 
     265         304 :                 m_aoStack.pop();
     266         304 :                 if (m_aoStack.empty())
     267         153 :                     break;
     268         151 :                 continue;
     269             :             }
     270        4127 :             auto &topContext = m_aoStack.top();
     271        4127 :             auto &sCurrentEntry = topContext.sEntries[topContext.nIdxInEntries];
     272        4127 :             if (sCurrentEntry.run_length == 0)
     273             :             {
     274             :                 // Arbitrary limit. 5 seems to be the maximum value supported
     275             :                 // by pmtiles.hpp::get_tile()
     276         215 :                 if (m_aoStack.size() == 5)
     277             :                 {
     278           0 :                     m_bEOF = true;
     279           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     280             :                              "Too many levels of nested directories");
     281           0 :                     break;
     282             :                 }
     283             : 
     284         430 :                 if (sHeader.leaf_dirs_offset >
     285         215 :                     std::numeric_limits<uint64_t>::max() - sCurrentEntry.offset)
     286             :                 {
     287           0 :                     m_bEOF = true;
     288           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     289             :                              "Invalid directory offset");
     290           0 :                     break;
     291             :                 }
     292         430 :                 const auto *posStr = m_poDS->ReadInternal(
     293         215 :                     sHeader.leaf_dirs_offset + sCurrentEntry.offset,
     294         215 :                     sCurrentEntry.length, "directory");
     295         215 :                 if (!posStr)
     296             :                 {
     297           0 :                     m_bEOF = true;
     298           0 :                     CPLError(
     299             :                         CE_Failure, CPLE_AppDefined,
     300             :                         "PMTILES: cannot read directory of size " CPL_FRMT_GUIB
     301             :                         " at "
     302             :                         "offset " CPL_FRMT_GUIB,
     303           0 :                         static_cast<GUIntBig>(sHeader.leaf_dirs_offset +
     304           0 :                                               sCurrentEntry.offset),
     305           0 :                         static_cast<GUIntBig>(sCurrentEntry.length));
     306           0 :                     break;
     307             :                 }
     308             : 
     309         215 :                 DirectoryContext sContext;
     310         215 :                 sContext.sEntries = pmtiles::deserialize_directory(*posStr);
     311         215 :                 if (sContext.sEntries.empty())
     312             :                 {
     313           0 :                     m_bEOF = true;
     314             :                     // In theory empty directories could exist, but for now
     315             :                     // do not allow this to be more robust against hostile files
     316             :                     // that could create many such empty directories
     317           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     318             :                              "Empty directory found");
     319           0 :                     break;
     320             :                 }
     321             : 
     322         317 :                 if (m_nLastTileId != INVALID_LAST_TILE_ID &&
     323         102 :                     sContext.sEntries[0].tile_id <= m_nLastTileId)
     324             :                 {
     325           0 :                     m_bEOF = true;
     326           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     327             :                              "Non increasing tile_id");
     328           0 :                     break;
     329             :                 }
     330             : 
     331         215 :                 if (m_nZoomLevel >= 0)
     332             :                 {
     333         215 :                     const int nMinEntryIdx = find_tile_idx_lesser_or_equal(
     334             :                         sContext.sEntries, m_nMinTileId);
     335         215 :                     if (nMinEntryIdx < 0)
     336             :                     {
     337           0 :                         if (AdvanceToNextTile())
     338           0 :                             continue;
     339           0 :                         break;
     340             :                     }
     341         215 :                     sContext.nIdxInEntries = nMinEntryIdx;
     342             :                 }
     343         215 :                 m_nLastTileId =
     344         215 :                     sContext.sEntries[sContext.nIdxInEntries].tile_id;
     345             : 
     346         215 :                 m_aoStack.emplace(std::move(sContext));
     347             : 
     348         215 :                 topContext.nIdxInEntries++;
     349             :             }
     350             :             else
     351             :             {
     352        3912 :                 if (topContext.nIdxInRunLength == sCurrentEntry.run_length)
     353             :                 {
     354        1152 :                     topContext.nIdxInEntries++;
     355        1152 :                     topContext.nIdxInRunLength = 0;
     356             :                 }
     357             :                 else
     358             :                 {
     359        2760 :                     const auto nIdxInRunLength = topContext.nIdxInRunLength;
     360        2760 :                     const uint64_t nTileId =
     361        2760 :                         sCurrentEntry.tile_id + nIdxInRunLength;
     362        2760 :                     m_nLastTileId = nTileId;
     363        2760 :                     const pmtiles::zxy zxy = pmtiles::tileid_to_zxy(nTileId);
     364             : 
     365             :                     // Sanity check to limit risk of iterating forever on
     366             :                     // broken run_length value
     367        2789 :                     if (nIdxInRunLength == 0 && sCurrentEntry.run_length > 1 &&
     368          29 :                         sCurrentEntry.run_length >
     369          29 :                             pmtiles::zxy_to_tileid(zxy.z + 1, 0, 0) - nTileId)
     370             :                     {
     371             :                         // Hit on /vsis3/us-west-2.opendata.source.coop/protomaps/openstreetmap/v4.pmtiles
     372           0 :                         sCurrentEntry.run_length = static_cast<uint32_t>(
     373           0 :                             pmtiles::zxy_to_tileid(zxy.z + 1, 0, 0) - nTileId);
     374             :                     }
     375             : 
     376        2760 :                     topContext.nIdxInRunLength++;
     377             : 
     378        2760 :                     if (m_nZoomLevel >= 0)
     379             :                     {
     380        2760 :                         if (nTileId < m_nMinTileId)
     381             :                         {
     382         843 :                             if (sCurrentEntry.run_length > 1)
     383             :                             {
     384          19 :                                 if (sCurrentEntry.tile_id +
     385          19 :                                         sCurrentEntry.run_length <=
     386          19 :                                     m_nMinTileId)
     387             :                                 {
     388           0 :                                     topContext.nIdxInRunLength =
     389           0 :                                         sCurrentEntry.run_length;
     390             :                                 }
     391             :                                 else
     392             :                                 {
     393          19 :                                     topContext.nIdxInRunLength =
     394             :                                         static_cast<uint32_t>(
     395          19 :                                             m_nMinTileId -
     396          19 :                                             sCurrentEntry.tile_id);
     397             :                                 }
     398          19 :                                 m_nLastTileId = sCurrentEntry.tile_id +
     399          19 :                                                 topContext.nIdxInRunLength - 1;
     400             :                             }
     401        1873 :                             continue;
     402             :                         }
     403        1917 :                         else if (nTileId > m_nMaxTileId)
     404             :                         {
     405         572 :                             if (AdvanceToNextTile())
     406         535 :                                 continue;
     407          37 :                             break;
     408             :                         }
     409             : 
     410        1345 :                         if (m_nMinX >= 0 &&
     411        1284 :                             zxy.x < static_cast<uint32_t>(m_nMinX))
     412         173 :                             continue;
     413        1172 :                         if (m_nMinY >= 0 &&
     414        1083 :                             zxy.y < static_cast<uint32_t>(m_nMinY))
     415          42 :                             continue;
     416        1130 :                         if (m_nMaxX >= 0 &&
     417        1069 :                             zxy.x > static_cast<uint32_t>(m_nMaxX))
     418         154 :                             continue;
     419         976 :                         if (m_nMaxY >= 0 &&
     420         887 :                             zxy.y > static_cast<uint32_t>(m_nMaxY))
     421         126 :                             continue;
     422             :                     }
     423             : 
     424        1700 :                     if (sHeader.tile_data_offset >
     425         850 :                         std::numeric_limits<uint64_t>::max() -
     426         850 :                             sCurrentEntry.offset)
     427             :                     {
     428           0 :                         m_bEOF = true;
     429           0 :                         CPLError(CE_Failure, CPLE_AppDefined,
     430             :                                  "Invalid tile offset");
     431           0 :                         break;
     432             :                     }
     433             : 
     434         850 :                     if (pnRunLength)
     435             :                     {
     436          18 :                         *pnRunLength =
     437          18 :                             sCurrentEntry.run_length - nIdxInRunLength;
     438             :                     }
     439             : 
     440             :                     // We must capture the result, before the below code
     441             :                     // that updates (m_nCurX, m_nCurY) and invalidates
     442             :                     // sCurrentEntry
     443             :                     const auto res = pmtiles::entry_zxy(
     444         850 :                         zxy.z, zxy.x, zxy.y,
     445         850 :                         sHeader.tile_data_offset + sCurrentEntry.offset,
     446         850 :                         sCurrentEntry.length);
     447             : 
     448         850 :                     AdvanceToNextTile();
     449             : 
     450         850 :                     return res;
     451             :                 }
     452             :             }
     453        3520 :         }
     454             :     }
     455           4 :     catch (const std::exception &e)
     456             :     {
     457           2 :         CPLError(CE_Failure, CPLE_AppDefined, "GetNextTile() failed with %s",
     458           2 :                  e.what());
     459             :     }
     460             : 
     461         192 :     m_bEOF = true;
     462         192 :     return pmtiles::entry_zxy(0, 0, 0, 0, 0);
     463             : }
     464             : 
     465             : /************************************************************************/
     466             : /*                             DumpTiles()                              */
     467             : /************************************************************************/
     468             : 
     469             : #ifdef DEBUG_PMTILES
     470             : void OGRPMTilesTileIterator::DumpTiles()
     471             : {
     472             :     int count = 0;
     473             :     while (true)
     474             :     {
     475             :         const auto sTile = GetNextTile();
     476             :         if (sTile.offset == 0)
     477             :             break;
     478             :         ++count;
     479             :         printf("%d -> %d %d %d %d %d\n",  // ok
     480             :                count, int(sTile.z), int(sTile.x), int(sTile.y),
     481             :                int(sTile.offset), int(sTile.length));
     482             :     }
     483             : }
     484             : #endif