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 : * Copyright (c) 2026, Even Rouault <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "ogr_pmtiles.h"
15 :
16 : #include "cpl_json.h"
17 :
18 : #include "mvtutils.h"
19 :
20 : #include <algorithm>
21 : #include <math.h>
22 :
23 : /************************************************************************/
24 : /* ~OGRPMTilesDataset() */
25 : /************************************************************************/
26 :
27 414 : OGRPMTilesDataset::~OGRPMTilesDataset()
28 : {
29 207 : if (!m_osMetadataFilename.empty())
30 185 : VSIUnlink(m_osMetadataFilename.c_str());
31 207 : if (!m_osTmpGPKGFilename.empty())
32 56 : VSIUnlink(m_osTmpGPKGFilename.c_str());
33 414 : }
34 :
35 : /************************************************************************/
36 : /* GetLayer() */
37 : /************************************************************************/
38 :
39 113 : const OGRLayer *OGRPMTilesDataset::GetLayer(int iLayer) const
40 :
41 : {
42 113 : if (iLayer < 0 || iLayer >= GetLayerCount())
43 10 : return nullptr;
44 103 : return m_apoLayers[iLayer].get();
45 : }
46 :
47 : /************************************************************************/
48 : /* LongLatToSphericalMercator() */
49 : /************************************************************************/
50 :
51 142 : static void LongLatToSphericalMercator(double *x, double *y)
52 : {
53 142 : double X = SPHERICAL_RADIUS * (*x) / 180 * M_PI;
54 142 : double Y = SPHERICAL_RADIUS * log(tan(M_PI / 4 + 0.5 * (*y) / 180 * M_PI));
55 142 : *x = X;
56 142 : *y = Y;
57 142 : }
58 :
59 : /************************************************************************/
60 : /* GetCompression() */
61 : /************************************************************************/
62 :
63 273 : /*static*/ const char *OGRPMTilesDataset::GetCompression(uint8_t nVal)
64 : {
65 273 : switch (nVal)
66 : {
67 6 : case pmtiles::COMPRESSION_UNKNOWN:
68 6 : return "unknown";
69 38 : case pmtiles::COMPRESSION_NONE:
70 38 : return "none";
71 229 : case pmtiles::COMPRESSION_GZIP:
72 229 : return "gzip";
73 0 : case pmtiles::COMPRESSION_BROTLI:
74 0 : return "brotli";
75 0 : case pmtiles::COMPRESSION_ZSTD:
76 0 : return "zstd";
77 0 : default:
78 0 : break;
79 : }
80 0 : return CPLSPrintf("invalid (%d)", nVal);
81 : }
82 :
83 : /************************************************************************/
84 : /* GetTileType() */
85 : /************************************************************************/
86 :
87 : /* static */
88 107 : const char *OGRPMTilesDataset::GetTileType(const pmtiles::headerv3 &sHeader)
89 : {
90 107 : switch (sHeader.tile_type)
91 : {
92 0 : case pmtiles::TILETYPE_UNKNOWN:
93 0 : return "unknown";
94 4 : case pmtiles::TILETYPE_MVT:
95 4 : return "MVT";
96 57 : case pmtiles::TILETYPE_PNG:
97 57 : return "PNG";
98 22 : case pmtiles::TILETYPE_JPEG:
99 22 : return "JPEG";
100 24 : case pmtiles::TILETYPE_WEBP:
101 24 : return "WEBP";
102 0 : case pmtiles::TILETYPE_AVIF:
103 0 : return "AVIF";
104 0 : case pmtiles::TILETYPE_MAPLIBRE_VECTOR_TILE:
105 0 : return "MAPLIBRE_VECTOR_TILE";
106 0 : default:
107 0 : break;
108 : }
109 0 : return CPLSPrintf("invalid (%d)", sHeader.tile_type);
110 : }
111 :
112 : /************************************************************************/
113 : /* Open() */
114 : /************************************************************************/
115 :
116 189 : bool OGRPMTilesDataset::Open(GDALOpenInfo *poOpenInfo)
117 : {
118 189 : if (!poOpenInfo->fpL || poOpenInfo->nHeaderBytes < PMTILES_HEADER_LENGTH)
119 2 : return false;
120 :
121 187 : SetDescription(poOpenInfo->pszFilename);
122 :
123 : // Borrow file handle
124 187 : m_poFileUniquePtr.reset(poOpenInfo->fpL);
125 187 : poOpenInfo->fpL = nullptr;
126 :
127 187 : m_poFile = m_poFileUniquePtr.get();
128 :
129 : // Deserizalize header
130 374 : std::string osHeader;
131 187 : osHeader.assign(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
132 187 : PMTILES_HEADER_LENGTH);
133 : try
134 : {
135 187 : m_sHeader = pmtiles::deserialize_header(osHeader);
136 : }
137 1 : catch (const std::exception &)
138 : {
139 1 : return false;
140 : }
141 :
142 : // Check tile type
143 186 : const bool bAcceptAnyTileType = CPLTestBool(CSLFetchNameValueDef(
144 186 : poOpenInfo->papszOpenOptions, "ACCEPT_ANY_TILE_TYPE", "NO"));
145 186 : if (bAcceptAnyTileType)
146 : {
147 : // do nothing. Internal use only by /vsipmtiles/
148 : }
149 79 : else if (m_sHeader.tile_type == pmtiles::TILETYPE_MVT)
150 : {
151 41 : if ((poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) == 0)
152 : {
153 0 : CPLError(CE_Failure, CPLE_AppDefined,
154 : "Tile type %s not handled in raster mode",
155 0 : GetTileType(m_sHeader));
156 0 : return false;
157 : }
158 : }
159 38 : else if (m_sHeader.tile_type == pmtiles::TILETYPE_PNG ||
160 15 : m_sHeader.tile_type == pmtiles::TILETYPE_JPEG ||
161 7 : m_sHeader.tile_type == pmtiles::TILETYPE_WEBP ||
162 0 : m_sHeader.tile_type == pmtiles::TILETYPE_AVIF)
163 : {
164 38 : if ((poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0)
165 : {
166 0 : CPLError(CE_Failure, CPLE_AppDefined,
167 : "Tile type %s not handled in vector mode",
168 0 : GetTileType(m_sHeader));
169 0 : return false;
170 : }
171 : }
172 : else
173 : {
174 0 : CPLError(CE_Failure, CPLE_AppDefined,
175 : "Tile type %s not handled by the driver",
176 0 : GetTileType(m_sHeader));
177 0 : return false;
178 : }
179 :
180 : // Check compression method for metadata and directories
181 186 : CPLDebugOnly("PMTiles", "internal_compression = %s",
182 : GetCompression(m_sHeader.internal_compression));
183 :
184 186 : if (m_sHeader.internal_compression == pmtiles::COMPRESSION_GZIP)
185 : {
186 186 : m_psInternalDecompressor = CPLGetDecompressor("gzip");
187 : }
188 0 : else if (m_sHeader.internal_compression == pmtiles::COMPRESSION_ZSTD)
189 : {
190 0 : m_psInternalDecompressor = CPLGetDecompressor("zstd");
191 0 : if (m_psInternalDecompressor == nullptr)
192 : {
193 0 : CPLError(CE_Failure, CPLE_AppDefined,
194 : "File %s requires ZSTD decompression, but not available "
195 : "in this GDAL build",
196 : poOpenInfo->pszFilename);
197 0 : return false;
198 : }
199 : }
200 0 : else if (m_sHeader.internal_compression != pmtiles::COMPRESSION_NONE)
201 : {
202 0 : CPLError(CE_Failure, CPLE_AppDefined,
203 : "Unhandled internal_compression = %s",
204 0 : GetCompression(m_sHeader.internal_compression));
205 0 : return false;
206 : }
207 :
208 : // Check compression for tile data
209 186 : if (!CPLTestBool(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
210 : "DECOMPRESS_TILES", "YES")))
211 : {
212 : // do nothing. Internal use only by /vsipmtiles/
213 : }
214 : else
215 : {
216 79 : CPLDebugOnly("PMTiles", "tile_compression = %s",
217 : GetCompression(m_sHeader.tile_compression));
218 :
219 79 : if (m_sHeader.tile_compression == pmtiles::COMPRESSION_UNKNOWN)
220 : {
221 : // Python pmtiles-convert generates this. The MVT driver can autodetect
222 : // uncompressed and GZip-compressed tiles automatically.
223 : }
224 75 : else if (m_sHeader.tile_compression == pmtiles::COMPRESSION_GZIP)
225 : {
226 37 : m_psTileDataDecompressor = CPLGetDecompressor("gzip");
227 : }
228 38 : else if (m_sHeader.tile_compression == pmtiles::COMPRESSION_ZSTD)
229 : {
230 0 : m_psTileDataDecompressor = CPLGetDecompressor("zstd");
231 0 : if (m_psTileDataDecompressor == nullptr)
232 : {
233 0 : CPLError(
234 : CE_Failure, CPLE_AppDefined,
235 : "File %s requires ZSTD decompression, but not available "
236 : "in this GDAL build",
237 : poOpenInfo->pszFilename);
238 0 : return false;
239 : }
240 : }
241 38 : else if (m_sHeader.tile_compression != pmtiles::COMPRESSION_NONE)
242 : {
243 0 : CPLError(CE_Failure, CPLE_AppDefined,
244 : "Unhandled tile_compression = %s",
245 0 : GetCompression(m_sHeader.tile_compression));
246 0 : return false;
247 : }
248 : }
249 :
250 : // Read metadata
251 : const auto *posMetadata =
252 186 : ReadInternal(m_sHeader.json_metadata_offset,
253 : m_sHeader.json_metadata_bytes, "metadata");
254 186 : if (!posMetadata)
255 1 : return false;
256 185 : CPLDebugOnly("PMTiles", "Metadata = %s", posMetadata->c_str());
257 185 : m_osMetadata = *posMetadata;
258 :
259 : m_osMetadataFilename =
260 185 : VSIMemGenerateHiddenFilename("pmtiles_metadata.json");
261 185 : VSIFCloseL(VSIFileFromMemBuffer(m_osMetadataFilename.c_str(),
262 185 : reinterpret_cast<GByte *>(&m_osMetadata[0]),
263 185 : m_osMetadata.size(), false));
264 :
265 370 : CPLJSONDocument oJsonDoc;
266 185 : if (!oJsonDoc.LoadMemory(m_osMetadata))
267 : {
268 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot parse metadata");
269 0 : return false;
270 : }
271 :
272 370 : auto oJsonRoot = oJsonDoc.GetRoot();
273 2015 : for (const auto &oChild : oJsonRoot.GetChildren())
274 : {
275 1830 : if (oChild.GetType() == CPLJSONObject::Type::String)
276 : {
277 1712 : if (oChild.GetName() == "json")
278 : {
279 : // Tippecanoe metadata includes a "json" item, which is a
280 : // serialized JSON object with vector_layers[] and layers[]
281 : // arrays we are interested in later.
282 : // so use "json" content as the new root
283 33 : if (!oJsonDoc.LoadMemory(oChild.ToString()))
284 : {
285 0 : CPLError(CE_Failure, CPLE_AppDefined,
286 : "Cannot parse 'json' metadata item");
287 0 : return false;
288 : }
289 33 : oJsonRoot = oJsonDoc.GetRoot();
290 : }
291 : // Tippecanoe generates a "strategies" member with serialized JSON
292 1679 : else if (oChild.GetName() != "strategies")
293 : {
294 1679 : SetMetadataItem(oChild.GetName().c_str(),
295 3358 : oChild.ToString().c_str());
296 : }
297 : }
298 : }
299 :
300 185 : m_nMinZoomLevel = m_sHeader.min_zoom;
301 185 : m_nMaxZoomLevel = m_sHeader.max_zoom;
302 185 : if (m_nMinZoomLevel > m_nMaxZoomLevel)
303 : {
304 1 : CPLError(CE_Failure, CPLE_AppDefined, "min_zoom(=%d) > max_zoom(=%d)",
305 : m_nMinZoomLevel, m_nMaxZoomLevel);
306 1 : return false;
307 : }
308 184 : if (m_nMinZoomLevel > 30)
309 : {
310 0 : CPLError(CE_Warning, CPLE_AppDefined, "Clamping min_zoom from %d to %d",
311 : m_nMinZoomLevel, 30);
312 0 : m_nMinZoomLevel = 30;
313 : }
314 184 : if (m_nMaxZoomLevel > 30)
315 : {
316 0 : CPLError(CE_Warning, CPLE_AppDefined, "Clamping max_zoom from %d to %d",
317 : m_nMaxZoomLevel, 30);
318 0 : m_nMaxZoomLevel = 30;
319 : }
320 :
321 184 : if (bAcceptAnyTileType)
322 106 : return true;
323 :
324 : const int nZoomLevel =
325 78 : atoi(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "ZOOM_LEVEL",
326 : CPLSPrintf("%d", m_nMaxZoomLevel)));
327 78 : if (nZoomLevel < m_nMinZoomLevel || nZoomLevel > m_nMaxZoomLevel)
328 : {
329 2 : CPLError(CE_Failure, CPLE_AppDefined,
330 : "Invalid zoom level. Should be in [%d,%d] range",
331 : m_nMinZoomLevel, m_nMaxZoomLevel);
332 2 : return false;
333 : }
334 76 : SetMetadataItem("ZOOM_LEVEL", CPLSPrintf("%d", nZoomLevel));
335 :
336 : m_osClipOpenOption =
337 76 : CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "CLIP", "");
338 :
339 76 : if (m_sHeader.tile_type == pmtiles::TILETYPE_MVT)
340 38 : return OpenVector(poOpenInfo, oJsonRoot, nZoomLevel);
341 : else
342 38 : return OpenRaster(nZoomLevel);
343 : }
344 :
345 : /************************************************************************/
346 : /* OpenVector() */
347 : /************************************************************************/
348 :
349 38 : bool OGRPMTilesDataset::OpenVector(const GDALOpenInfo *poOpenInfo,
350 : const CPLJSONObject &oJsonRoot,
351 : int nZoomLevel)
352 : {
353 : // If using the pmtiles go utility, vector_layers and tilestats are
354 : // moved from Tippecanoe's json metadata item to the root element.
355 114 : CPLJSONArray oVectorLayers = oJsonRoot.GetArray("vector_layers");
356 38 : if (oVectorLayers.Size() == 0)
357 : {
358 5 : CPLError(CE_Failure, CPLE_AppDefined,
359 : "Missing vector_layers[] metadata");
360 5 : return false;
361 : }
362 :
363 66 : CPLJSONArray oTileStatLayers = oJsonRoot.GetArray("tilestats/layers");
364 :
365 66 : const bool bZoomLevelFromSpatialFilter = CPLFetchBool(
366 33 : poOpenInfo->papszOpenOptions, "ZOOM_LEVEL_AUTO",
367 33 : CPLTestBool(CPLGetConfigOption("MVT_ZOOM_LEVEL_AUTO", "NO")));
368 : const bool bJsonField =
369 33 : CPLFetchBool(poOpenInfo->papszOpenOptions, "JSON_FIELD", false);
370 :
371 33 : double dfMinX = m_sHeader.min_lon_e7 / 10e6;
372 33 : double dfMinY = m_sHeader.min_lat_e7 / 10e6;
373 33 : double dfMaxX = m_sHeader.max_lon_e7 / 10e6;
374 33 : double dfMaxY = m_sHeader.max_lat_e7 / 10e6;
375 33 : LongLatToSphericalMercator(&dfMinX, &dfMinY);
376 33 : LongLatToSphericalMercator(&dfMaxX, &dfMaxY);
377 :
378 73 : for (int i = 0; i < oVectorLayers.Size(); i++)
379 : {
380 120 : CPLJSONObject oId = oVectorLayers[i].GetObj("id");
381 40 : if (oId.IsValid() && oId.GetType() == CPLJSONObject::Type::String)
382 : {
383 40 : OGRwkbGeometryType eGeomType = wkbUnknown;
384 40 : if (oTileStatLayers.IsValid())
385 : {
386 38 : eGeomType = OGRMVTFindGeomTypeFromTileStat(
387 76 : oTileStatLayers, oId.ToString().c_str());
388 : }
389 40 : if (eGeomType == wkbUnknown)
390 : {
391 2 : eGeomType = OGRPMTilesVectorLayer::GuessGeometryType(
392 4 : this, oId.ToString().c_str(), nZoomLevel);
393 : }
394 :
395 120 : CPLJSONObject oFields = oVectorLayers[i].GetObj("fields");
396 : CPLJSONArray oAttributesFromTileStats =
397 : OGRMVTFindAttributesFromTileStat(oTileStatLayers,
398 80 : oId.ToString().c_str());
399 :
400 40 : m_apoLayers.push_back(std::make_unique<OGRPMTilesVectorLayer>(
401 80 : this, oId.ToString().c_str(), oFields, oAttributesFromTileStats,
402 : bJsonField, dfMinX, dfMinY, dfMaxX, dfMaxY, eGeomType,
403 : nZoomLevel, bZoomLevelFromSpatialFilter));
404 : }
405 : }
406 :
407 33 : return true;
408 : }
409 :
410 : /************************************************************************/
411 : /* OpenRaster() */
412 : /************************************************************************/
413 :
414 38 : bool OGRPMTilesDataset::OpenRaster(int nZoomLevel)
415 : {
416 38 : m_nZoomLevel = nZoomLevel;
417 :
418 38 : m_oSRS.importFromEPSG(3857);
419 :
420 : // Open a tile to get its dimension in pixel (to compute the resolution)
421 76 : OGRPMTilesTileIterator oIter(this, nZoomLevel);
422 38 : auto sTile = oIter.GetNextTile();
423 38 : if (sTile.offset == 0)
424 : {
425 0 : CPLError(CE_Failure, CPLE_AppDefined,
426 : "File does not contain any tile at zoom level %d", nZoomLevel);
427 0 : return false;
428 : }
429 :
430 38 : const auto *posStr = ReadTileData(sTile.offset, sTile.length);
431 38 : if (!posStr)
432 : {
433 : // Error message emitted by ReadTileData()
434 0 : return false;
435 : }
436 38 : const auto &osTileData = *posStr;
437 :
438 38 : auto poDM = GetGDALDriverManager();
439 38 : if (!poDM->GetDriverByName(GetTileType(m_sHeader)))
440 : {
441 0 : CPLError(CE_Failure, CPLE_AppDefined,
442 : "%s driver needed for inner working of PMTiles driver is not "
443 : "available",
444 0 : GetTileType(m_sHeader));
445 0 : return false;
446 : }
447 38 : m_poGPKGDriver = poDM->GetDriverByName("GPKG");
448 38 : if (!m_poGPKGDriver)
449 : {
450 0 : CPLError(CE_Failure, CPLE_AppDefined,
451 : "GeoPackage driver needed for inner working of PMTiles driver "
452 : "is not available");
453 0 : return false;
454 : }
455 38 : if (!poDM->GetDriverByName("GTI"))
456 : {
457 0 : CPLError(CE_Failure, CPLE_AppDefined,
458 : "GTI driver needed for inner working of PMTiles driver is not "
459 : "available");
460 0 : return false;
461 : }
462 : const std::string osTmpFilename = VSIMemGenerateHiddenFilename(
463 76 : CPLSPrintf("pmtiles_%u_%u_tile", sTile.x, sTile.y));
464 38 : VSIFCloseL(VSIFileFromMemBuffer(
465 : osTmpFilename.c_str(),
466 38 : reinterpret_cast<GByte *>(const_cast<char *>(osTileData.data())),
467 38 : osTileData.size(), false));
468 :
469 38 : const char *const apszAllowedDrivers[] = {GetTileType(m_sHeader), nullptr};
470 : auto poTileDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
471 : osTmpFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
472 76 : apszAllowedDrivers));
473 38 : const int nTileSize = poTileDS ? poTileDS->GetRasterXSize() : 0;
474 38 : if (nTileSize == 0)
475 0 : return false;
476 :
477 38 : const double dfRes = 2 * MAX_GM / (1 << nZoomLevel) / nTileSize;
478 38 : constexpr double EPSILON = 1e-2;
479 :
480 : // Compute the raster georeferenced extent from the bounding box in the
481 : // PMTiles metadata
482 38 : double dfMinX = m_sHeader.min_lon_e7 / 10e6;
483 38 : double dfMinY = m_sHeader.min_lat_e7 / 10e6;
484 38 : double dfMaxX = m_sHeader.max_lon_e7 / 10e6;
485 38 : double dfMaxY = m_sHeader.max_lat_e7 / 10e6;
486 38 : LongLatToSphericalMercator(&dfMinX, &dfMinY);
487 38 : LongLatToSphericalMercator(&dfMaxX, &dfMaxY);
488 :
489 : // Align with the resolution at the zoom level of interest, to avoid
490 : // resampling due to sub-pixel shifts
491 38 : dfMinX = std::floor(dfMinX / dfRes + EPSILON) * dfRes;
492 38 : dfMinY = std::floor(dfMinY / dfRes + EPSILON) * dfRes;
493 38 : dfMaxX = std::ceil(dfMaxX / dfRes - EPSILON) * dfRes;
494 38 : dfMaxY = std::ceil(dfMaxY / dfRes - EPSILON) * dfRes;
495 :
496 : // Compute the geotransform
497 38 : m_gt.xorig = dfMinX;
498 38 : m_gt.xscale = dfRes;
499 38 : m_gt.yorig = dfMaxY;
500 38 : m_gt.yscale = -dfRes;
501 :
502 : // Compute the raster dimension
503 38 : const double dfXSize = (dfMaxX - dfMinX) / dfRes;
504 38 : const double dfYSize = (dfMaxY - dfMinY) / dfRes;
505 38 : if (dfXSize > INT_MAX || dfYSize > INT_MAX)
506 : {
507 0 : CPLError(CE_Failure, CPLE_AppDefined,
508 : "Too large zoom level %d compared to dataset extent",
509 : m_nZoomLevel);
510 0 : return false;
511 : }
512 38 : nRasterXSize = std::max(1, static_cast<int>(dfXSize + EPSILON));
513 38 : nRasterYSize = std::max(1, static_cast<int>(dfYSize + EPSILON));
514 :
515 : // Create bands
516 38 : const int nTargetBands = poTileDS->GetRasterBand(1)->GetColorTable()
517 38 : ? 3
518 38 : : poTileDS->GetRasterCount();
519 180 : for (int i = 0; i < nTargetBands; ++i)
520 142 : SetBand(i + 1, std::make_unique<GDALPMTilesRasterBand>(this, i + 1,
521 : nTileSize));
522 :
523 38 : m_osTmpGPKGFilename = VSIMemGenerateHiddenFilename("pmtiles.gti.gpkg");
524 :
525 38 : if (nTargetBands > 1)
526 38 : SetMetadataItem("INTERLEAVING", "PIXEL", "IMAGE_STRUCTURE");
527 :
528 : // Create overview datasets
529 56 : for (int iOvr = 0; iOvr < m_nZoomLevel - m_nMinZoomLevel; ++iOvr)
530 : {
531 36 : auto poOvrDS = std::make_unique<OGRPMTilesDataset>();
532 18 : poOvrDS->SetDescription(GetDescription());
533 18 : poOvrDS->m_psInternalDecompressor = m_psInternalDecompressor;
534 18 : poOvrDS->m_psTileDataDecompressor = m_psTileDataDecompressor;
535 18 : poOvrDS->m_poFile = m_poFile;
536 18 : poOvrDS->m_nZoomLevel = m_nZoomLevel - 1 - iOvr;
537 18 : poOvrDS->m_osTmpGPKGFilename = m_osTmpGPKGFilename;
538 18 : poOvrDS->m_poGPKGDriver = m_poGPKGDriver;
539 18 : poOvrDS->m_sHeader = m_sHeader;
540 18 : poOvrDS->m_oSRS = m_oSRS;
541 18 : poOvrDS->m_gt = m_gt;
542 18 : poOvrDS->m_gt.xscale *= (1 << (iOvr + 1));
543 18 : poOvrDS->m_gt.yscale *= (1 << (iOvr + 1));
544 18 : poOvrDS->nRasterXSize = std::max(
545 36 : 1,
546 18 : static_cast<int>((dfMaxX - dfMinX) / poOvrDS->m_gt.xscale + 0.5));
547 18 : poOvrDS->nRasterYSize = std::max(
548 36 : 1,
549 18 : static_cast<int>((dfMaxY - dfMinY) / -poOvrDS->m_gt.yscale + 0.5));
550 18 : poOvrDS->m_gt.xscale = (dfMaxX - dfMinX) / poOvrDS->nRasterXSize;
551 18 : poOvrDS->m_gt.yscale = -(dfMaxY - dfMinY) / poOvrDS->nRasterYSize;
552 :
553 86 : for (int i = 0; i < nTargetBands; ++i)
554 136 : poOvrDS->SetBand(i + 1, std::make_unique<GDALPMTilesRasterBand>(
555 136 : poOvrDS.get(), i + 1, nTileSize));
556 :
557 18 : if (nTargetBands > 1)
558 18 : poOvrDS->SetMetadataItem("INTERLEAVING", "PIXEL",
559 18 : "IMAGE_STRUCTURE");
560 :
561 18 : m_apoOverviews.push_back(std::move(poOvrDS));
562 : }
563 :
564 38 : return true;
565 : }
566 :
567 : /************************************************************************/
568 : /* IRasterIO() */
569 : /************************************************************************/
570 :
571 28 : CPLErr OGRPMTilesDataset::IRasterIO(
572 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
573 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
574 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
575 : GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
576 : {
577 :
578 28 : if (nBufXSize < nXSize && nBufYSize < nYSize && AreOverviewsEnabled())
579 : {
580 1 : int bTried = FALSE;
581 1 : const CPLErr eErr = TryOverviewRasterIO(
582 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
583 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
584 : nBandSpace, psExtraArg, &bTried);
585 1 : if (bTried)
586 1 : return eErr;
587 : }
588 :
589 : // Compute the georeferenced coordinates of the window of interest
590 : // defined by (nXOff, nYOff, nXSize, nYSize) (or the floating point
591 : // coordinates)
592 27 : const double dfXOff =
593 27 : psExtraArg->bFloatingPointWindowValidity ? psExtraArg->dfXOff : nXOff;
594 27 : const double dfYOff =
595 27 : psExtraArg->bFloatingPointWindowValidity ? psExtraArg->dfYOff : nYOff;
596 27 : const double dfXSize =
597 27 : psExtraArg->bFloatingPointWindowValidity ? psExtraArg->dfXSize : nXSize;
598 27 : const double dfYSize =
599 27 : psExtraArg->bFloatingPointWindowValidity ? psExtraArg->dfYSize : nYSize;
600 27 : const double dfMinX = m_gt.xorig + dfXOff * m_gt.xscale;
601 27 : const double dfMaxY = m_gt.yorig + dfYOff * m_gt.yscale;
602 27 : const double dfMaxX = m_gt.xorig + (dfXOff + dfXSize) * m_gt.xscale;
603 27 : const double dfMinY = m_gt.yorig + (dfYOff + dfYSize) * m_gt.yscale;
604 27 : const double dfTileDim = 2 * MAX_GM / (1 << m_nZoomLevel);
605 :
606 : // Compute the minimum and maximum tile indices covering the window
607 : // of interest
608 27 : constexpr double EPSILON = 1e-5;
609 : const int nTileMinX = std::max(
610 27 : 0, static_cast<int>(floor((dfMinX + MAX_GM) / dfTileDim + EPSILON)));
611 : // PMTiles uses a Y=MAX_GM as the y=0 tile
612 : const int nTileMinY = std::max(
613 27 : 0, static_cast<int>(floor((MAX_GM - dfMaxY) / dfTileDim + EPSILON)));
614 : const int nTileMaxX = std::min(
615 54 : static_cast<int>(floor((dfMaxX + MAX_GM) / dfTileDim + EPSILON)),
616 27 : (1 << m_nZoomLevel) - 1);
617 : const int nTileMaxY = std::min(
618 54 : static_cast<int>(floor((MAX_GM - dfMinY) / dfTileDim + EPSILON)),
619 27 : (1 << m_nZoomLevel) - 1);
620 :
621 : // Create a GTI GeoPackage file with the tiles that intersect the window
622 : // of interest.
623 27 : VSIUnlink(m_osTmpGPKGFilename.c_str());
624 27 : auto poGPKGDS = std::unique_ptr<GDALDataset>(m_poGPKGDriver->Create(
625 54 : m_osTmpGPKGFilename.c_str(), 0, 0, 0, GDT_Unknown, nullptr));
626 27 : auto poLayer = poGPKGDS ? poGPKGDS->CreateLayer("index") : nullptr;
627 27 : if (!poLayer)
628 : {
629 0 : return CE_Failure;
630 : }
631 54 : OGRFieldDefn oFieldDefn("location", OFTString);
632 27 : CPL_IGNORE_RET_VAL(poLayer->CreateField(&oFieldDefn));
633 :
634 : // Add layer-level metadata items recognized by the GTI driver.
635 27 : const double dfMosaicMinX = -MAX_GM + nTileMinX * dfTileDim;
636 27 : const double dfMosaicMinY = MAX_GM - (nTileMaxY + 1) * dfTileDim;
637 27 : const double dfMosaicMaxX = -MAX_GM + (nTileMaxX + 1) * dfTileDim;
638 27 : const double dfMosaicMaxY = MAX_GM - nTileMinY * dfTileDim;
639 27 : poLayer->SetMetadataItem("RESX", CPLSPrintf("%.17g", m_gt.xscale));
640 27 : poLayer->SetMetadataItem("RESY", CPLSPrintf("%.17g", -m_gt.yscale));
641 27 : poLayer->SetMetadataItem("MINX", CPLSPrintf("%.17g", dfMosaicMinX));
642 27 : poLayer->SetMetadataItem("MINY", CPLSPrintf("%.17g", dfMosaicMinY));
643 27 : poLayer->SetMetadataItem("MAXX", CPLSPrintf("%.17g", dfMosaicMaxX));
644 27 : poLayer->SetMetadataItem("MAXY", CPLSPrintf("%.17g", dfMosaicMaxY));
645 27 : poLayer->SetMetadataItem("DATA_TYPE", "UInt8");
646 27 : poLayer->SetMetadataItem("BAND_COUNT", CPLSPrintf("%d", nBands));
647 27 : if (nBands <= 4)
648 : {
649 27 : poLayer->SetMetadataItem("COLOR_INTERPRETATION",
650 27 : nBands == 1 ? "gray"
651 49 : : nBands == 2 ? "gray,alpha"
652 22 : : nBands == 3 ? "red,green,blue"
653 27 : : "red,green,blue,alpha");
654 : }
655 :
656 27 : poGPKGDS->StartTransaction();
657 : OGRPMTilesTileIterator oTileIter(this, m_nZoomLevel, nTileMinX, nTileMinY,
658 54 : nTileMaxX, nTileMaxY);
659 : while (true)
660 : {
661 54 : const auto sTile = oTileIter.GetNextTile();
662 54 : if (sTile.offset == 0)
663 : {
664 27 : break;
665 : }
666 :
667 54 : OGRFeature oFeature(poLayer->GetLayerDefn());
668 54 : oFeature.SetField(0, CPLSPrintf("/vsipmtiles/%s/%d/%d/%d%s",
669 27 : GetDescription(), m_nZoomLevel, sTile.x,
670 27 : sTile.y,
671 : VSIPMTilesGetTileExtension(this)));
672 54 : auto poLR = std::make_unique<OGRLinearRing>();
673 27 : const double dfTileMinX = -MAX_GM + sTile.x * dfTileDim;
674 27 : const double dfTileMaxX = -MAX_GM + (sTile.x + 1) * dfTileDim;
675 27 : const double dfTileMaxY = MAX_GM - sTile.y * dfTileDim;
676 27 : const double dfTileMinY = MAX_GM - (sTile.y + 1) * dfTileDim;
677 27 : poLR->addPoint(dfTileMinX, dfTileMinY);
678 27 : poLR->addPoint(dfTileMinX, dfTileMaxY);
679 27 : poLR->addPoint(dfTileMaxX, dfTileMaxY);
680 27 : poLR->addPoint(dfTileMaxX, dfTileMinY);
681 27 : poLR->addPoint(dfTileMinX, dfTileMinY);
682 27 : auto poPoly = std::make_unique<OGRPolygon>();
683 27 : poPoly->addRing(std::move(poLR));
684 27 : oFeature.SetGeometry(std::move(poPoly));
685 27 : CPL_IGNORE_RET_VAL(poLayer->CreateFeature(&oFeature));
686 27 : }
687 27 : poGPKGDS->CommitTransaction();
688 27 : poGPKGDS.reset();
689 :
690 : // Open the GTI GeoPackage file has a raster
691 27 : const char *const apszAllowedDriversGTI[] = {"GTI", nullptr};
692 54 : CPLStringList aosOpenOptionsGTI;
693 : aosOpenOptionsGTI.SetNameValue("ALLOWED_RASTER_DRIVERS",
694 27 : GetTileType(m_sHeader));
695 : auto poGTIDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
696 : m_osTmpGPKGFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
697 54 : apszAllowedDriversGTI, aosOpenOptionsGTI.List()));
698 27 : if (!poGTIDS)
699 : {
700 0 : return CE_Failure;
701 : }
702 :
703 : // Compute the window of interest relative to the origin of the GTI dataset
704 : GDALRasterIOExtraArg sExtraArg;
705 27 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
706 27 : sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
707 27 : sExtraArg.bFloatingPointWindowValidity =
708 27 : psExtraArg->bFloatingPointWindowValidity;
709 27 : sExtraArg.dfXOff = (dfMinX - dfMosaicMinX) / m_gt.xscale;
710 27 : sExtraArg.dfYOff = (dfMosaicMaxY - dfMaxY) / -m_gt.yscale;
711 27 : sExtraArg.dfXSize = psExtraArg->dfXSize;
712 27 : sExtraArg.dfYSize = psExtraArg->dfYSize;
713 :
714 27 : const int nGTIXOff = static_cast<int>(sExtraArg.dfXOff + EPSILON);
715 27 : const int nGTIYOff = static_cast<int>(sExtraArg.dfYOff + EPSILON);
716 :
717 : // Redirect the pixel request to the GTI dataset
718 27 : return poGTIDS->RasterIO(GF_Read, nGTIXOff, nGTIYOff, nXSize, nYSize, pData,
719 : nBufXSize, nBufYSize, eBufType, nBandCount,
720 : panBandMap, nPixelSpace, nLineSpace, nBandSpace,
721 27 : &sExtraArg);
722 : }
723 :
724 : /************************************************************************/
725 : /* Read() */
726 : /************************************************************************/
727 :
728 1100 : const std::string *OGRPMTilesDataset::Read(const CPLCompressor *psDecompressor,
729 : uint64_t nOffset, uint64_t nSize,
730 : const char *pszDataType)
731 : {
732 1100 : if (nSize > 10 * 1024 * 1024)
733 : {
734 0 : CPLError(CE_Failure, CPLE_AppDefined,
735 : "Too large amount of %s to read: " CPL_FRMT_GUIB
736 : " bytes at offset " CPL_FRMT_GUIB,
737 : pszDataType, static_cast<GUIntBig>(nSize),
738 : static_cast<GUIntBig>(nOffset));
739 0 : return nullptr;
740 : }
741 1100 : m_osBuffer.resize(static_cast<size_t>(nSize));
742 2200 : if (m_poFile->Seek(nOffset, SEEK_SET) != 0 ||
743 1100 : m_poFile->Read(&m_osBuffer[0], m_osBuffer.size(), 1) != 1)
744 : {
745 1 : CPLError(CE_Failure, CPLE_AppDefined,
746 : "Cannot read %s of length %u at offset " CPL_FRMT_GUIB,
747 : pszDataType, unsigned(nSize), static_cast<GUIntBig>(nOffset));
748 1 : return nullptr;
749 : }
750 :
751 1099 : if (psDecompressor)
752 : {
753 917 : m_osDecompressedBuffer.resize(32 + 16 * m_osBuffer.size());
754 917 : for (int iTry = 0; iTry < 2; ++iTry)
755 : {
756 917 : void *pOutputData = &m_osDecompressedBuffer[0];
757 917 : size_t nOutputSize = m_osDecompressedBuffer.size();
758 917 : if (!psDecompressor->pfnFunc(m_osBuffer.data(), m_osBuffer.size(),
759 : &pOutputData, &nOutputSize, nullptr,
760 917 : psDecompressor->user_data))
761 : {
762 0 : if (iTry == 0)
763 : {
764 0 : pOutputData = nullptr;
765 0 : nOutputSize = 0;
766 0 : if (psDecompressor->pfnFunc(
767 0 : m_osBuffer.data(), m_osBuffer.size(), &pOutputData,
768 0 : &nOutputSize, nullptr, psDecompressor->user_data))
769 : {
770 0 : CPLDebug("PMTiles",
771 : "Buffer of size %u uncompresses to %u bytes",
772 : unsigned(nSize), unsigned(nOutputSize));
773 0 : m_osDecompressedBuffer.resize(nOutputSize);
774 0 : continue;
775 : }
776 : }
777 :
778 0 : CPLError(CE_Failure, CPLE_AppDefined,
779 : "Cannot decompress %s of length %u at "
780 : "offset " CPL_FRMT_GUIB,
781 : pszDataType, unsigned(nSize),
782 : static_cast<GUIntBig>(nOffset));
783 0 : return nullptr;
784 : }
785 917 : m_osDecompressedBuffer.resize(nOutputSize);
786 917 : break;
787 : }
788 917 : return &m_osDecompressedBuffer;
789 : }
790 : else
791 : {
792 182 : return &m_osBuffer;
793 : }
794 : }
795 :
796 : /************************************************************************/
797 : /* ReadInternal() */
798 : /************************************************************************/
799 :
800 778 : const std::string *OGRPMTilesDataset::ReadInternal(uint64_t nOffset,
801 : uint64_t nSize,
802 : const char *pszDataType)
803 : {
804 778 : return Read(m_psInternalDecompressor, nOffset, nSize, pszDataType);
805 : }
806 :
807 : /************************************************************************/
808 : /* ReadTileData() */
809 : /************************************************************************/
810 :
811 322 : const std::string *OGRPMTilesDataset::ReadTileData(uint64_t nOffset,
812 : uint64_t nSize)
813 : {
814 322 : return Read(m_psTileDataDecompressor, nOffset, nSize, "tile data");
815 : }
|
