Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL Rasterlite driver
4 : * Purpose: Implement GDAL Rasterlite support using OGR SQLite driver
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : **********************************************************************
8 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_string.h"
14 : #include "gdal_frmts.h"
15 : #include "ogr_api.h"
16 : #include "ogr_srs_api.h"
17 :
18 : #include "rasterlitedataset.h"
19 : #include "rasterlitedrivercore.h"
20 :
21 : #include <algorithm>
22 : #include <limits>
23 :
24 : /************************************************************************/
25 : /* RasterliteOpenSQLiteDB() */
26 : /************************************************************************/
27 :
28 84 : GDALDatasetH RasterliteOpenSQLiteDB(const char *pszFilename, GDALAccess eAccess)
29 : {
30 84 : const char *const apszAllowedDrivers[] = {"SQLITE", nullptr};
31 84 : return GDALOpenEx(pszFilename,
32 : GDAL_OF_VECTOR |
33 : ((eAccess == GA_Update) ? GDAL_OF_UPDATE : 0),
34 168 : apszAllowedDrivers, nullptr, nullptr);
35 : }
36 :
37 : /************************************************************************/
38 : /* RasterliteGetPixelSizeCond() */
39 : /************************************************************************/
40 :
41 104 : CPLString RasterliteGetPixelSizeCond(double dfPixelXSize, double dfPixelYSize,
42 : const char *pszTablePrefixWithDot)
43 : {
44 104 : CPLString osCond;
45 : osCond.Printf("((%spixel_x_size >= %s AND %spixel_x_size <= %s) AND "
46 : "(%spixel_y_size >= %s AND %spixel_y_size <= %s))",
47 : pszTablePrefixWithDot,
48 208 : CPLString().FormatC(dfPixelXSize - 1e-15, "%.15f").c_str(),
49 : pszTablePrefixWithDot,
50 208 : CPLString().FormatC(dfPixelXSize + 1e-15, "%.15f").c_str(),
51 : pszTablePrefixWithDot,
52 208 : CPLString().FormatC(dfPixelYSize - 1e-15, "%.15f").c_str(),
53 : pszTablePrefixWithDot,
54 520 : CPLString().FormatC(dfPixelYSize + 1e-15, "%.15f").c_str());
55 104 : return osCond;
56 : }
57 :
58 : /************************************************************************/
59 : /* RasterliteGetSpatialFilterCond() */
60 : /************************************************************************/
61 :
62 45 : CPLString RasterliteGetSpatialFilterCond(double minx, double miny, double maxx,
63 : double maxy)
64 : {
65 45 : CPLString osCond;
66 : osCond.Printf("(xmin < %s AND xmax > %s AND ymin < %s AND ymax > %s)",
67 90 : CPLString().FormatC(maxx, "%.15f").c_str(),
68 90 : CPLString().FormatC(minx, "%.15f").c_str(),
69 90 : CPLString().FormatC(maxy, "%.15f").c_str(),
70 225 : CPLString().FormatC(miny, "%.15f").c_str());
71 45 : return osCond;
72 : }
73 :
74 : /************************************************************************/
75 : /* RasterliteBand() */
76 : /************************************************************************/
77 :
78 76 : RasterliteBand::RasterliteBand(RasterliteDataset *poDSIn, int nBandIn,
79 : GDALDataType eDataTypeIn, int nBlockXSizeIn,
80 76 : int nBlockYSizeIn)
81 : {
82 76 : poDS = poDSIn;
83 76 : nBand = nBandIn;
84 76 : eDataType = eDataTypeIn;
85 76 : nBlockXSize = nBlockXSizeIn;
86 76 : nBlockYSize = nBlockYSizeIn;
87 76 : }
88 :
89 : /************************************************************************/
90 : /* IReadBlock() */
91 : /************************************************************************/
92 :
93 : // #define RASTERLITE_DEBUG
94 :
95 18 : CPLErr RasterliteBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
96 : {
97 18 : RasterliteDataset *poGDS = reinterpret_cast<RasterliteDataset *>(poDS);
98 :
99 18 : double minx = poGDS->adfGeoTransform[0] +
100 18 : nBlockXOff * nBlockXSize * poGDS->adfGeoTransform[1];
101 18 : double maxx = poGDS->adfGeoTransform[0] +
102 18 : (nBlockXOff + 1) * nBlockXSize * poGDS->adfGeoTransform[1];
103 18 : double maxy = poGDS->adfGeoTransform[3] +
104 18 : nBlockYOff * nBlockYSize * poGDS->adfGeoTransform[5];
105 18 : double miny = poGDS->adfGeoTransform[3] +
106 18 : (nBlockYOff + 1) * nBlockYSize * poGDS->adfGeoTransform[5];
107 18 : int nDataTypeSize = GDALGetDataTypeSize(eDataType) / 8;
108 :
109 : #ifdef RASTERLITE_DEBUG
110 : if (nBand == 1)
111 : {
112 : printf("nBlockXOff = %d, nBlockYOff = %d, nBlockXSize = %d, " /*ok*/
113 : "nBlockYSize = %d\n"
114 : "minx=%.15f miny=%.15f maxx=%.15f maxy=%.15f\n",
115 : nBlockXOff, nBlockYOff, nBlockXSize, nBlockYSize, minx, miny,
116 : maxx, maxy);
117 : }
118 : #endif
119 :
120 36 : CPLString osSQL;
121 : osSQL.Printf("SELECT m.geometry, r.raster, m.id, m.width, m.height FROM "
122 : "\"%s_metadata\" AS m, "
123 : "\"%s_rasters\" AS r WHERE m.rowid IN "
124 : "(SELECT pkid FROM \"idx_%s_metadata_geometry\" "
125 : "WHERE %s) AND %s AND r.id = m.id",
126 : poGDS->osTableName.c_str(), poGDS->osTableName.c_str(),
127 : poGDS->osTableName.c_str(),
128 36 : RasterliteGetSpatialFilterCond(minx, miny, maxx, maxy).c_str(),
129 18 : RasterliteGetPixelSizeCond(poGDS->adfGeoTransform[1],
130 18 : -poGDS->adfGeoTransform[5], "m.")
131 36 : .c_str());
132 :
133 : OGRLayerH hSQLLyr =
134 18 : GDALDatasetExecuteSQL(poGDS->hDS, osSQL.c_str(), nullptr, nullptr);
135 18 : if (hSQLLyr == nullptr)
136 : {
137 0 : memset(pImage, 0,
138 0 : static_cast<size_t>(nBlockXSize) * nBlockYSize * nDataTypeSize);
139 0 : return CE_None;
140 : }
141 :
142 : const CPLString osMemFileName(
143 36 : VSIMemGenerateHiddenFilename("rasterlite_tile"));
144 :
145 : #ifdef RASTERLITE_DEBUG
146 : if (nBand == 1)
147 : {
148 : printf("nTiles = %d\n", /*ok*/
149 : static_cast<int>(OGR_L_GetFeatureCount(hSQLLyr, TRUE)));
150 : }
151 : #endif
152 :
153 18 : bool bHasFoundTile = false;
154 18 : bool bHasMemsetTile = false;
155 :
156 : OGRFeatureH hFeat;
157 18 : CPLErr eErr = CE_None;
158 36 : while ((hFeat = OGR_L_GetNextFeature(hSQLLyr)) != nullptr &&
159 : eErr == CE_None)
160 : {
161 18 : OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat);
162 18 : if (hGeom == nullptr)
163 : {
164 0 : CPLError(CE_Failure, CPLE_AppDefined, "null geometry found");
165 0 : OGR_F_Destroy(hFeat);
166 0 : GDALDatasetReleaseResultSet(poGDS->hDS, hSQLLyr);
167 0 : memset(pImage, 0,
168 0 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
169 0 : nDataTypeSize);
170 0 : return CE_Failure;
171 : }
172 :
173 18 : OGREnvelope oEnvelope;
174 18 : OGR_G_GetEnvelope(hGeom, &oEnvelope);
175 :
176 18 : const int nTileId = OGR_F_GetFieldAsInteger(hFeat, 1);
177 18 : if (poGDS->m_nLastBadTileId == nTileId)
178 : {
179 0 : OGR_F_Destroy(hFeat);
180 0 : continue;
181 : }
182 18 : const int nTileXSize = OGR_F_GetFieldAsInteger(hFeat, 2);
183 18 : const int nTileYSize = OGR_F_GetFieldAsInteger(hFeat, 3);
184 18 : constexpr int MAX_INT_DIV_2 = std::numeric_limits<int>::max() / 2;
185 18 : if (nTileXSize <= 0 || nTileXSize >= MAX_INT_DIV_2 || nTileYSize <= 0 ||
186 : nTileYSize >= MAX_INT_DIV_2)
187 : {
188 0 : CPLError(CE_Failure, CPLE_AppDefined, "invalid tile size");
189 0 : OGR_F_Destroy(hFeat);
190 0 : GDALDatasetReleaseResultSet(poGDS->hDS, hSQLLyr);
191 0 : memset(pImage, 0,
192 0 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
193 0 : nDataTypeSize);
194 0 : return CE_Failure;
195 : }
196 :
197 18 : const double dfDstXOff =
198 18 : (oEnvelope.MinX - minx) / poGDS->adfGeoTransform[1];
199 18 : const double dfDstYOff =
200 18 : (maxy - oEnvelope.MaxY) / (-poGDS->adfGeoTransform[5]);
201 18 : constexpr int MIN_INT_DIV_2 = std::numeric_limits<int>::min() / 2;
202 18 : if (!(dfDstXOff >= MIN_INT_DIV_2 && dfDstXOff <= MAX_INT_DIV_2) ||
203 18 : !(dfDstYOff >= MIN_INT_DIV_2 && dfDstYOff <= MAX_INT_DIV_2))
204 : {
205 0 : CPLError(CE_Failure, CPLE_AppDefined, "invalid geometry");
206 0 : OGR_F_Destroy(hFeat);
207 0 : GDALDatasetReleaseResultSet(poGDS->hDS, hSQLLyr);
208 0 : memset(pImage, 0,
209 0 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
210 0 : nDataTypeSize);
211 0 : return CE_Failure;
212 : }
213 18 : int nDstXOff = static_cast<int>(dfDstXOff + 0.5);
214 18 : int nDstYOff = static_cast<int>(dfDstYOff + 0.5);
215 :
216 18 : int nReqXSize = nTileXSize;
217 18 : int nReqYSize = nTileYSize;
218 :
219 : int nSrcXOff, nSrcYOff;
220 :
221 18 : if (nDstXOff >= 0)
222 : {
223 17 : nSrcXOff = 0;
224 : }
225 : else
226 : {
227 1 : nSrcXOff = -nDstXOff;
228 1 : nReqXSize += nDstXOff;
229 1 : nDstXOff = 0;
230 : }
231 :
232 18 : if (nDstYOff >= 0)
233 : {
234 17 : nSrcYOff = 0;
235 : }
236 : else
237 : {
238 1 : nSrcYOff = -nDstYOff;
239 1 : nReqYSize += nDstYOff;
240 1 : nDstYOff = 0;
241 : }
242 :
243 18 : if (nDstXOff + nReqXSize > nBlockXSize)
244 0 : nReqXSize = nBlockXSize - nDstXOff;
245 :
246 18 : if (nDstYOff + nReqYSize > nBlockYSize)
247 0 : nReqYSize = nBlockYSize - nDstYOff;
248 :
249 : #ifdef RASTERLITE_DEBUG
250 : if (nBand == 1)
251 : {
252 : printf(/*ok*/
253 : "id = %d, minx=%.15f miny=%.15f maxx=%.15f maxy=%.15f\n"
254 : "nDstXOff = %d, nDstYOff = %d, nSrcXOff = %d, nSrcYOff = "
255 : "%d, "
256 : "nReqXSize=%d, nReqYSize=%d\n",
257 : nTileId, oEnvelope.MinX, oEnvelope.MinY, oEnvelope.MaxX,
258 : oEnvelope.MaxY, nDstXOff, nDstYOff, nSrcXOff, nSrcYOff,
259 : nReqXSize, nReqYSize);
260 : }
261 : #endif
262 :
263 18 : if (nReqXSize > 0 && nReqYSize > 0 && nSrcXOff < nTileXSize &&
264 : nSrcYOff < nTileYSize)
265 : {
266 :
267 : #ifdef RASTERLITE_DEBUG
268 : if (nBand == 1)
269 : {
270 : printf("id = %d, selected !\n", nTileId); /*ok*/
271 : }
272 : #endif
273 18 : int nDataSize = 0;
274 18 : GByte *pabyData = OGR_F_GetFieldAsBinary(hFeat, 0, &nDataSize);
275 :
276 18 : VSILFILE *fp = VSIFileFromMemBuffer(osMemFileName.c_str(), pabyData,
277 : nDataSize, FALSE);
278 18 : VSIFCloseL(fp);
279 :
280 18 : GDALDatasetH hDSTile = GDALOpenEx(osMemFileName.c_str(),
281 : GDAL_OF_RASTER | GDAL_OF_INTERNAL,
282 : nullptr, nullptr, nullptr);
283 18 : int nTileBands = 0;
284 18 : if (hDSTile && (nTileBands = GDALGetRasterCount(hDSTile)) == 0)
285 : {
286 0 : GDALClose(hDSTile);
287 0 : hDSTile = nullptr;
288 : }
289 18 : if (hDSTile == nullptr)
290 : {
291 0 : poGDS->m_nLastBadTileId = nTileId;
292 0 : CPLError(CE_Failure, CPLE_AppDefined, "Can't open tile %d",
293 : nTileId);
294 : }
295 :
296 18 : int nReqBand = 1;
297 18 : if (nTileBands == poGDS->nBands)
298 17 : nReqBand = nBand;
299 1 : else if (eDataType == GDT_Byte && nTileBands == 1 &&
300 1 : poGDS->nBands == 3)
301 1 : nReqBand = 1;
302 : else
303 : {
304 0 : poGDS->m_nLastBadTileId = nTileId;
305 0 : GDALClose(hDSTile);
306 0 : hDSTile = nullptr;
307 : }
308 :
309 18 : if (hDSTile)
310 : {
311 36 : if (GDALGetRasterXSize(hDSTile) != nTileXSize ||
312 18 : GDALGetRasterYSize(hDSTile) != nTileYSize)
313 : {
314 0 : CPLError(CE_Failure, CPLE_AppDefined,
315 : "Invalid dimensions for tile %d", nTileId);
316 0 : poGDS->m_nLastBadTileId = nTileId;
317 0 : GDALClose(hDSTile);
318 0 : hDSTile = nullptr;
319 : }
320 : }
321 :
322 18 : if (hDSTile)
323 : {
324 18 : bHasFoundTile = true;
325 :
326 18 : bool bHasJustMemsetTileBand1 = false;
327 :
328 : /* If the source tile doesn't fit the entire block size, then */
329 : /* we memset 0 before */
330 18 : if (!(nDstXOff == 0 && nDstYOff == 0 &&
331 18 : nReqXSize == nBlockXSize && nReqYSize == nBlockYSize) &&
332 6 : !bHasMemsetTile)
333 : {
334 6 : memset(pImage, 0,
335 6 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
336 6 : nDataTypeSize);
337 6 : bHasMemsetTile = true;
338 6 : bHasJustMemsetTileBand1 = true;
339 : }
340 :
341 : GDALColorTable *poTileCT = reinterpret_cast<GDALColorTable *>(
342 18 : GDALGetRasterColorTable(GDALGetRasterBand(hDSTile, 1)));
343 18 : unsigned char *pabyTranslationTable = nullptr;
344 18 : if (poGDS->nBands == 1 && poGDS->poCT != nullptr &&
345 : poTileCT != nullptr)
346 : {
347 : pabyTranslationTable = reinterpret_cast<GDALRasterBand *>(
348 1 : GDALGetRasterBand(hDSTile, 1))
349 1 : ->GetIndexColorTranslationTo(
350 : this, nullptr, nullptr);
351 : }
352 :
353 : /* --------------------------------------------------------------------
354 : */
355 : /* Read tile data */
356 : /* --------------------------------------------------------------------
357 : */
358 18 : eErr = GDALRasterIO(
359 : GDALGetRasterBand(hDSTile, nReqBand), GF_Read, nSrcXOff,
360 : nSrcYOff, nReqXSize, nReqYSize,
361 : reinterpret_cast<char *>(pImage) +
362 18 : (nDstXOff + nDstYOff * nBlockXSize) * nDataTypeSize,
363 : nReqXSize, nReqYSize, eDataType, nDataTypeSize,
364 18 : nBlockXSize * nDataTypeSize);
365 :
366 18 : if (eDataType == GDT_Byte && pabyTranslationTable)
367 : {
368 : /* --------------------------------------------------------------------
369 : */
370 : /* Convert from tile CT to band CT */
371 : /* --------------------------------------------------------------------
372 : */
373 0 : for (int j = nDstYOff; j < nDstYOff + nReqYSize; j++)
374 : {
375 0 : for (int i = nDstXOff; i < nDstXOff + nReqXSize; i++)
376 : {
377 0 : GByte *pPixel = reinterpret_cast<GByte *>(pImage) +
378 0 : i + j * nBlockXSize;
379 0 : *pPixel = pabyTranslationTable[*pPixel];
380 : }
381 : }
382 0 : CPLFree(pabyTranslationTable);
383 0 : pabyTranslationTable = nullptr;
384 : }
385 18 : else if (eDataType == GDT_Byte && nTileBands == 1 &&
386 15 : poGDS->nBands == 3 && poTileCT != nullptr)
387 : {
388 : /* --------------------------------------------------------------------
389 : */
390 : /* Expand from PCT to RGB */
391 : /* --------------------------------------------------------------------
392 : */
393 : GByte abyCT[256];
394 : const int nEntries =
395 1 : std::min(256, poTileCT->GetColorEntryCount());
396 257 : for (int i = 0; i < nEntries; i++)
397 : {
398 : const GDALColorEntry *psEntry =
399 256 : poTileCT->GetColorEntry(i);
400 256 : if (nBand == 1)
401 256 : abyCT[i] = static_cast<GByte>(psEntry->c1);
402 0 : else if (nBand == 2)
403 0 : abyCT[i] = static_cast<GByte>(psEntry->c2);
404 : else
405 0 : abyCT[i] = static_cast<GByte>(psEntry->c3);
406 : }
407 1 : for (int i = nEntries; i < 256; i++)
408 0 : abyCT[i] = 0;
409 :
410 170 : for (int j = nDstYOff; j < nDstYOff + nReqYSize; j++)
411 : {
412 57291 : for (int i = nDstXOff; i < nDstXOff + nReqXSize; i++)
413 : {
414 57122 : GByte *pPixel = reinterpret_cast<GByte *>(pImage) +
415 57122 : i + j * nBlockXSize;
416 57122 : *pPixel = abyCT[*pPixel];
417 : }
418 : }
419 : }
420 :
421 : /* --------------------------------------------------------------------
422 : */
423 : /* Put in the block cache the data for this block into
424 : * other bands */
425 : /* while the underlying dataset is opened */
426 : /* --------------------------------------------------------------------
427 : */
428 18 : if (nBand == 1 && poGDS->nBands > 1)
429 : {
430 12 : for (int iOtherBand = 2;
431 12 : iOtherBand <= poGDS->nBands && eErr == CE_None;
432 : iOtherBand++)
433 : {
434 : GDALRasterBlock *poBlock =
435 8 : poGDS->GetRasterBand(iOtherBand)
436 16 : ->GetLockedBlockRef(nBlockXOff, nBlockYOff,
437 8 : TRUE);
438 8 : if (poBlock == nullptr)
439 0 : break;
440 :
441 : GByte *pabySrcBlock =
442 8 : reinterpret_cast<GByte *>(poBlock->GetDataRef());
443 8 : if (pabySrcBlock == nullptr)
444 : {
445 0 : poBlock->DropLock();
446 0 : break;
447 : }
448 :
449 8 : if (nTileBands == 1)
450 2 : nReqBand = 1;
451 : else
452 6 : nReqBand = iOtherBand;
453 :
454 8 : if (bHasJustMemsetTileBand1)
455 2 : memset(pabySrcBlock, 0,
456 2 : static_cast<size_t>(nBlockXSize) *
457 2 : nBlockYSize * nDataTypeSize);
458 :
459 : /* --------------------------------------------------------------------
460 : */
461 : /* Read tile data */
462 : /* --------------------------------------------------------------------
463 : */
464 8 : eErr = GDALRasterIO(
465 : GDALGetRasterBand(hDSTile, nReqBand), GF_Read,
466 : nSrcXOff, nSrcYOff, nReqXSize, nReqYSize,
467 : reinterpret_cast<char *>(pabySrcBlock) +
468 8 : (nDstXOff + nDstYOff * nBlockXSize) *
469 : nDataTypeSize,
470 : nReqXSize, nReqYSize, eDataType, nDataTypeSize,
471 8 : nBlockXSize * nDataTypeSize);
472 :
473 8 : if (eDataType == GDT_Byte && nTileBands == 1 &&
474 2 : poGDS->nBands == 3 && poTileCT != nullptr)
475 : {
476 : /* --------------------------------------------------------------------
477 : */
478 : /* Expand from PCT to RGB */
479 : /* --------------------------------------------------------------------
480 : */
481 :
482 : GByte abyCT[256];
483 : const int nEntries =
484 2 : std::min(256, poTileCT->GetColorEntryCount());
485 514 : for (int i = 0; i < nEntries; i++)
486 : {
487 : const GDALColorEntry *psEntry =
488 512 : poTileCT->GetColorEntry(i);
489 512 : if (iOtherBand == 2)
490 256 : abyCT[i] = static_cast<GByte>(psEntry->c2);
491 : else
492 256 : abyCT[i] = static_cast<GByte>(psEntry->c3);
493 : }
494 2 : for (int i = nEntries; i < 256; i++)
495 0 : abyCT[i] = 0;
496 :
497 340 : for (int j = nDstYOff; j < nDstYOff + nReqYSize;
498 : j++)
499 : {
500 114582 : for (int i = nDstXOff; i < nDstXOff + nReqXSize;
501 : i++)
502 : {
503 114244 : GByte *pPixel = reinterpret_cast<GByte *>(
504 : pabySrcBlock) +
505 114244 : i + j * nBlockXSize;
506 114244 : *pPixel = abyCT[*pPixel];
507 : }
508 : }
509 : }
510 :
511 8 : poBlock->DropLock();
512 : }
513 : }
514 18 : GDALClose(hDSTile);
515 : }
516 :
517 18 : VSIUnlink(osMemFileName.c_str());
518 : }
519 : else
520 : {
521 : #ifdef RASTERLITE_DEBUG
522 : if (nBand == 1)
523 : {
524 : printf("id = %d, NOT selected !\n", nTileId); /*ok*/
525 : }
526 : #endif
527 : }
528 18 : OGR_F_Destroy(hFeat);
529 : }
530 :
531 18 : VSIUnlink(osMemFileName.c_str());
532 18 : VSIUnlink((osMemFileName + ".aux.xml").c_str());
533 :
534 18 : if (!bHasFoundTile)
535 : {
536 0 : memset(pImage, 0,
537 0 : static_cast<size_t>(nBlockXSize) * nBlockYSize * nDataTypeSize);
538 : }
539 :
540 18 : GDALDatasetReleaseResultSet(poGDS->hDS, hSQLLyr);
541 :
542 : #ifdef RASTERLITE_DEBUG
543 : if (nBand == 1)
544 : printf("\n"); /*ok*/
545 : #endif
546 :
547 18 : return eErr;
548 : }
549 :
550 : /************************************************************************/
551 : /* GetOverviewCount() */
552 : /************************************************************************/
553 :
554 9 : int RasterliteBand::GetOverviewCount()
555 : {
556 9 : RasterliteDataset *poGDS = reinterpret_cast<RasterliteDataset *>(poDS);
557 :
558 9 : if (poGDS->nLimitOvrCount >= 0)
559 5 : return poGDS->nLimitOvrCount;
560 4 : else if (poGDS->nResolutions > 1)
561 1 : return poGDS->nResolutions - 1;
562 : else
563 3 : return GDALPamRasterBand::GetOverviewCount();
564 : }
565 :
566 : /************************************************************************/
567 : /* GetOverview() */
568 : /************************************************************************/
569 :
570 13 : GDALRasterBand *RasterliteBand::GetOverview(int nLevel)
571 : {
572 13 : RasterliteDataset *poGDS = reinterpret_cast<RasterliteDataset *>(poDS);
573 :
574 13 : if (poGDS->nLimitOvrCount >= 0)
575 : {
576 4 : if (nLevel < 0 || nLevel >= poGDS->nLimitOvrCount)
577 0 : return nullptr;
578 : }
579 :
580 13 : if (poGDS->nResolutions == 1)
581 0 : return GDALPamRasterBand::GetOverview(nLevel);
582 :
583 13 : if (nLevel < 0 || nLevel >= poGDS->nResolutions - 1)
584 0 : return nullptr;
585 :
586 13 : GDALDataset *poOvrDS = poGDS->papoOverviews[nLevel];
587 13 : if (poOvrDS)
588 13 : return poOvrDS->GetRasterBand(nBand);
589 :
590 0 : return nullptr;
591 : }
592 :
593 : /************************************************************************/
594 : /* GetColorInterpretation() */
595 : /************************************************************************/
596 :
597 1 : GDALColorInterp RasterliteBand::GetColorInterpretation()
598 : {
599 1 : RasterliteDataset *poGDS = reinterpret_cast<RasterliteDataset *>(poDS);
600 1 : if (poGDS->nBands == 1)
601 : {
602 1 : if (poGDS->poCT != nullptr)
603 1 : return GCI_PaletteIndex;
604 :
605 0 : return GCI_GrayIndex;
606 : }
607 0 : else if (poGDS->nBands == 3)
608 : {
609 0 : if (nBand == 1)
610 0 : return GCI_RedBand;
611 0 : else if (nBand == 2)
612 0 : return GCI_GreenBand;
613 0 : else if (nBand == 3)
614 0 : return GCI_BlueBand;
615 : }
616 :
617 0 : return GCI_Undefined;
618 : }
619 :
620 : /************************************************************************/
621 : /* GetColorTable() */
622 : /************************************************************************/
623 :
624 3 : GDALColorTable *RasterliteBand::GetColorTable()
625 : {
626 3 : RasterliteDataset *poGDS = reinterpret_cast<RasterliteDataset *>(poDS);
627 3 : if (poGDS->nBands == 1)
628 2 : return poGDS->poCT;
629 :
630 1 : return nullptr;
631 : }
632 :
633 : /************************************************************************/
634 : /* RasterliteDataset() */
635 : /************************************************************************/
636 :
637 72 : RasterliteDataset::RasterliteDataset()
638 : : bMustFree(FALSE), poMainDS(nullptr), nLevel(0), papszMetadata(nullptr),
639 144 : papszImageStructure(CSLAddString(nullptr, "INTERLEAVE=PIXEL")),
640 : papszSubDatasets(nullptr), nResolutions(0), padfXResolutions(nullptr),
641 : padfYResolutions(nullptr), papoOverviews(nullptr), nLimitOvrCount(-1),
642 : bValidGeoTransform(FALSE), poCT(nullptr), bCheckForExistingOverview(TRUE),
643 72 : hDS(nullptr)
644 : {
645 72 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
646 72 : memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
647 72 : }
648 :
649 : /************************************************************************/
650 : /* RasterliteDataset() */
651 : /************************************************************************/
652 :
653 15 : RasterliteDataset::RasterliteDataset(RasterliteDataset *poMainDSIn,
654 15 : int nLevelIn)
655 : : bMustFree(FALSE), poMainDS(poMainDSIn), nLevel(nLevelIn),
656 15 : papszMetadata(poMainDSIn->papszMetadata),
657 15 : papszImageStructure(poMainDSIn->papszImageStructure),
658 15 : papszSubDatasets(poMainDS->papszSubDatasets),
659 15 : nResolutions(poMainDSIn->nResolutions - nLevelIn),
660 15 : padfXResolutions(poMainDSIn->padfXResolutions + nLevelIn),
661 15 : padfYResolutions(poMainDSIn->padfYResolutions + nLevelIn),
662 15 : papoOverviews(poMainDSIn->papoOverviews + nLevelIn), nLimitOvrCount(-1),
663 15 : bValidGeoTransform(TRUE), m_oSRS(poMainDSIn->m_oSRS),
664 15 : poCT(poMainDSIn->poCT), osTableName(poMainDSIn->osTableName),
665 15 : osFileName(poMainDSIn->osFileName), bCheckForExistingOverview(TRUE),
666 : // TODO: osOvrFileName?
667 15 : hDS(poMainDSIn->hDS)
668 : {
669 15 : nRasterXSize = static_cast<int>(
670 15 : poMainDS->nRasterXSize *
671 15 : (poMainDS->padfXResolutions[0] / padfXResolutions[0]) +
672 : 0.5);
673 15 : nRasterYSize = static_cast<int>(
674 15 : poMainDS->nRasterYSize *
675 15 : (poMainDS->padfYResolutions[0] / padfYResolutions[0]) +
676 : 0.5);
677 :
678 15 : memcpy(adfGeoTransform, poMainDS->adfGeoTransform, 6 * sizeof(double));
679 15 : adfGeoTransform[1] = padfXResolutions[0];
680 15 : adfGeoTransform[5] = -padfYResolutions[0];
681 15 : }
682 :
683 : /************************************************************************/
684 : /* ~RasterliteDataset() */
685 : /************************************************************************/
686 :
687 174 : RasterliteDataset::~RasterliteDataset()
688 : {
689 87 : RasterliteDataset::CloseDependentDatasets();
690 174 : }
691 :
692 : /************************************************************************/
693 : /* CloseDependentDatasets() */
694 : /************************************************************************/
695 :
696 110 : int RasterliteDataset::CloseDependentDatasets()
697 : {
698 110 : int bRet = GDALPamDataset::CloseDependentDatasets();
699 :
700 110 : if (poMainDS == nullptr && !bMustFree)
701 : {
702 90 : CSLDestroy(papszMetadata);
703 90 : papszMetadata = nullptr;
704 90 : CSLDestroy(papszSubDatasets);
705 90 : papszSubDatasets = nullptr;
706 90 : CSLDestroy(papszImageStructure);
707 90 : papszImageStructure = nullptr;
708 :
709 90 : if (papoOverviews)
710 : {
711 18 : for (int i = 1; i < nResolutions; i++)
712 : {
713 11 : if (papoOverviews[i - 1] != nullptr &&
714 10 : papoOverviews[i - 1]->bMustFree)
715 : {
716 0 : papoOverviews[i - 1]->poMainDS = nullptr;
717 : }
718 11 : delete papoOverviews[i - 1];
719 : }
720 7 : CPLFree(papoOverviews);
721 7 : papoOverviews = nullptr;
722 7 : nResolutions = 0;
723 7 : bRet = TRUE;
724 : }
725 :
726 90 : if (hDS != nullptr)
727 42 : GDALClose(hDS);
728 90 : hDS = nullptr;
729 :
730 90 : CPLFree(padfXResolutions);
731 90 : CPLFree(padfYResolutions);
732 90 : padfXResolutions = nullptr;
733 90 : padfYResolutions = nullptr;
734 :
735 90 : delete poCT;
736 90 : poCT = nullptr;
737 : }
738 20 : else if (poMainDS != nullptr && bMustFree)
739 : {
740 1 : poMainDS->papoOverviews[nLevel - 1] = nullptr;
741 1 : delete poMainDS;
742 1 : poMainDS = nullptr;
743 1 : bRet = TRUE;
744 : }
745 :
746 110 : return bRet;
747 : }
748 :
749 : /************************************************************************/
750 : /* AddSubDataset() */
751 : /************************************************************************/
752 :
753 30 : void RasterliteDataset::AddSubDataset(const char *pszDSName)
754 : {
755 : char szName[80];
756 30 : const int nCount = CSLCount(papszSubDatasets) / 2;
757 :
758 30 : snprintf(szName, sizeof(szName), "SUBDATASET_%d_NAME", nCount + 1);
759 30 : papszSubDatasets = CSLSetNameValue(papszSubDatasets, szName, pszDSName);
760 :
761 30 : snprintf(szName, sizeof(szName), "SUBDATASET_%d_DESC", nCount + 1);
762 30 : papszSubDatasets = CSLSetNameValue(papszSubDatasets, szName, pszDSName);
763 30 : }
764 :
765 : /************************************************************************/
766 : /* GetMetadataDomainList() */
767 : /************************************************************************/
768 :
769 0 : char **RasterliteDataset::GetMetadataDomainList()
770 : {
771 0 : return BuildMetadataDomainList(GDALPamDataset::GetMetadataDomainList(),
772 : TRUE, "SUBDATASETS", "IMAGE_STRUCTURE",
773 0 : nullptr);
774 : }
775 :
776 : /************************************************************************/
777 : /* GetMetadata() */
778 : /************************************************************************/
779 :
780 15 : char **RasterliteDataset::GetMetadata(const char *pszDomain)
781 :
782 : {
783 15 : if (pszDomain != nullptr && EQUAL(pszDomain, "SUBDATASETS"))
784 0 : return papszSubDatasets;
785 :
786 15 : if (CSLCount(papszSubDatasets) < 2 && pszDomain != nullptr &&
787 0 : EQUAL(pszDomain, "IMAGE_STRUCTURE"))
788 0 : return papszImageStructure;
789 :
790 15 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
791 15 : return papszMetadata;
792 :
793 0 : return GDALPamDataset::GetMetadata(pszDomain);
794 : }
795 :
796 : /************************************************************************/
797 : /* GetMetadataItem() */
798 : /************************************************************************/
799 :
800 17 : const char *RasterliteDataset::GetMetadataItem(const char *pszName,
801 : const char *pszDomain)
802 : {
803 17 : if (pszDomain != nullptr && EQUAL(pszDomain, "OVERVIEWS"))
804 : {
805 2 : if (nResolutions > 1 || CSLCount(papszSubDatasets) > 2)
806 0 : return nullptr;
807 :
808 2 : osOvrFileName.Printf("%s_%s", osFileName.c_str(), osTableName.c_str());
809 4 : if (bCheckForExistingOverview == FALSE ||
810 2 : CPLCheckForFile(const_cast<char *>(osOvrFileName.c_str()), nullptr))
811 0 : return osOvrFileName.c_str();
812 :
813 2 : return nullptr;
814 : }
815 :
816 15 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
817 : }
818 :
819 : /************************************************************************/
820 : /* GetGeoTransform() */
821 : /************************************************************************/
822 :
823 2 : CPLErr RasterliteDataset::GetGeoTransform(double *padfGeoTransform)
824 : {
825 2 : if (bValidGeoTransform)
826 : {
827 2 : memcpy(padfGeoTransform, adfGeoTransform, 6 * sizeof(double));
828 2 : return CE_None;
829 : }
830 :
831 0 : return CE_Failure;
832 : }
833 :
834 : /************************************************************************/
835 : /* GetSpatialRef() */
836 : /************************************************************************/
837 :
838 2 : const OGRSpatialReference *RasterliteDataset::GetSpatialRef() const
839 : {
840 2 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
841 : }
842 :
843 : /************************************************************************/
844 : /* GetFileList() */
845 : /************************************************************************/
846 :
847 0 : char **RasterliteDataset::GetFileList()
848 : {
849 0 : char **papszFileList = CSLAddString(nullptr, osFileName);
850 0 : return papszFileList;
851 : }
852 :
853 : /************************************************************************/
854 : /* GetBlockParams() */
855 : /************************************************************************/
856 :
857 48 : int RasterliteDataset::GetBlockParams(OGRLayerH hRasterLyr, int nLevelIn,
858 : int *pnBands, GDALDataType *peDataType,
859 : int *pnBlockXSize, int *pnBlockYSize)
860 : {
861 96 : CPLString osSQL;
862 : osSQL.Printf("SELECT m.geometry, r.raster, m.id "
863 : "FROM \"%s_metadata\" AS m, \"%s_rasters\" AS r "
864 : "WHERE %s AND r.id = m.id",
865 : osTableName.c_str(), osTableName.c_str(),
866 48 : RasterliteGetPixelSizeCond(padfXResolutions[nLevelIn],
867 48 : padfYResolutions[nLevelIn], "m.")
868 48 : .c_str());
869 :
870 : OGRLayerH hSQLLyr =
871 48 : GDALDatasetExecuteSQL(hDS, osSQL.c_str(), nullptr, nullptr);
872 48 : if (hSQLLyr == nullptr)
873 : {
874 0 : return FALSE;
875 : }
876 :
877 48 : OGRFeatureH hFeat = OGR_L_GetNextFeature(hRasterLyr);
878 48 : if (hFeat == nullptr)
879 : {
880 0 : GDALDatasetReleaseResultSet(hDS, hSQLLyr);
881 0 : return FALSE;
882 : }
883 :
884 : int nDataSize;
885 48 : GByte *pabyData = OGR_F_GetFieldAsBinary(hFeat, 0, &nDataSize);
886 :
887 48 : if (nDataSize > 32 &&
888 48 : STARTS_WITH_CI(reinterpret_cast<const char *>(pabyData),
889 : "StartWaveletsImage$$"))
890 : {
891 0 : CPLError(CE_Failure, CPLE_NotSupported,
892 : "Rasterlite driver no longer support WAVELET compressed "
893 : "images");
894 0 : OGR_F_Destroy(hFeat);
895 0 : GDALDatasetReleaseResultSet(hDS, hSQLLyr);
896 0 : return FALSE;
897 : }
898 :
899 : const CPLString osMemFileName(
900 48 : VSIMemGenerateHiddenFilename("rasterlite_tile"));
901 : VSILFILE *fp =
902 48 : VSIFileFromMemBuffer(osMemFileName.c_str(), pabyData, nDataSize, FALSE);
903 48 : VSIFCloseL(fp);
904 :
905 48 : GDALDatasetH hDSTile = GDALOpen(osMemFileName.c_str(), GA_ReadOnly);
906 48 : if (hDSTile)
907 : {
908 48 : *pnBands = GDALGetRasterCount(hDSTile);
909 48 : if (*pnBands == 0)
910 : {
911 0 : GDALClose(hDSTile);
912 0 : hDSTile = nullptr;
913 : }
914 : }
915 : else
916 : {
917 0 : CPLError(CE_Failure, CPLE_AppDefined, "Can't open tile %d",
918 : OGR_F_GetFieldAsInteger(hFeat, 1));
919 : }
920 :
921 48 : if (hDSTile)
922 : {
923 48 : *peDataType = GDALGetRasterDataType(GDALGetRasterBand(hDSTile, 1));
924 :
925 74 : for (int iBand = 2; iBand <= *pnBands; iBand++)
926 : {
927 52 : if (*peDataType !=
928 26 : GDALGetRasterDataType(GDALGetRasterBand(hDSTile, 1)))
929 : {
930 0 : CPLError(CE_Failure, CPLE_NotSupported,
931 : "Band types must be identical");
932 0 : GDALClose(hDSTile);
933 0 : hDSTile = nullptr;
934 0 : goto end;
935 : }
936 : }
937 :
938 48 : *pnBlockXSize = GDALGetRasterXSize(hDSTile);
939 48 : *pnBlockYSize = GDALGetRasterYSize(hDSTile);
940 48 : if (CSLFindName(papszImageStructure, "COMPRESSION") == -1)
941 : {
942 : const char *pszCompression =
943 45 : GDALGetMetadataItem(hDSTile, "COMPRESSION", "IMAGE_STRUCTURE");
944 45 : if (pszCompression != nullptr && EQUAL(pszCompression, "JPEG"))
945 5 : papszImageStructure =
946 5 : CSLAddString(papszImageStructure, "COMPRESSION=JPEG");
947 : }
948 :
949 48 : if (CSLFindName(papszMetadata, "TILE_FORMAT") == -1)
950 : {
951 33 : papszMetadata = CSLSetNameValue(
952 : papszMetadata, "TILE_FORMAT",
953 : GDALGetDriverShortName(GDALGetDatasetDriver(hDSTile)));
954 : }
955 :
956 48 : if (*pnBands == 1 && this->poCT == nullptr)
957 : {
958 : GDALColorTable *l_poCT = reinterpret_cast<GDALColorTable *>(
959 36 : GDALGetRasterColorTable(GDALGetRasterBand(hDSTile, 1)));
960 36 : if (l_poCT)
961 2 : this->poCT = l_poCT->Clone();
962 : }
963 :
964 48 : GDALClose(hDSTile);
965 : }
966 0 : end:
967 48 : VSIUnlink(osMemFileName.c_str());
968 48 : VSIUnlink((osMemFileName + ".aux.xml").c_str());
969 :
970 48 : OGR_F_Destroy(hFeat);
971 :
972 48 : GDALDatasetReleaseResultSet(hDS, hSQLLyr);
973 :
974 48 : return hDSTile != nullptr;
975 : }
976 :
977 : /************************************************************************/
978 : /* Open() */
979 : /************************************************************************/
980 :
981 54 : GDALDataset *RasterliteDataset::Open(GDALOpenInfo *poOpenInfo)
982 : {
983 54 : if (RasterliteDriverIdentify(poOpenInfo) == FALSE)
984 0 : return nullptr;
985 :
986 108 : CPLString osFileName;
987 108 : CPLString osTableName;
988 54 : int nLevel = 0;
989 54 : double minx = 0.0;
990 54 : double miny = 0.0;
991 54 : double maxx = 0.0;
992 54 : double maxy = 0.0;
993 54 : int bMinXSet = FALSE;
994 54 : int bMinYSet = FALSE;
995 54 : int bMaxXSet = FALSE;
996 54 : int bMaxYSet = FALSE;
997 54 : int nReqBands = 0;
998 :
999 : /* -------------------------------------------------------------------- */
1000 : /* Parse "file name" */
1001 : /* -------------------------------------------------------------------- */
1002 : #ifdef ENABLE_SQL_SQLITE_FORMAT
1003 54 : if (poOpenInfo->pabyHeader &&
1004 37 : STARTS_WITH((const char *)poOpenInfo->pabyHeader, "-- SQL RASTERLITE"))
1005 : {
1006 1 : osFileName = poOpenInfo->pszFilename;
1007 : }
1008 : else
1009 : #endif
1010 53 : if (poOpenInfo->nHeaderBytes >= 1024 && poOpenInfo->pabyHeader &&
1011 36 : STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader,
1012 : "SQLite Format 3"))
1013 : {
1014 36 : osFileName = poOpenInfo->pszFilename;
1015 : }
1016 : else
1017 : {
1018 34 : char **papszTokens = CSLTokenizeStringComplex(
1019 17 : poOpenInfo->pszFilename + 11, ",", FALSE, FALSE);
1020 17 : int nTokens = CSLCount(papszTokens);
1021 17 : if (nTokens == 0)
1022 : {
1023 0 : CSLDestroy(papszTokens);
1024 0 : return nullptr;
1025 : }
1026 :
1027 17 : osFileName = papszTokens[0];
1028 :
1029 38 : for (int i = 1; i < nTokens; i++)
1030 : {
1031 21 : if (STARTS_WITH_CI(papszTokens[i], "table="))
1032 15 : osTableName = papszTokens[i] + 6;
1033 6 : else if (STARTS_WITH_CI(papszTokens[i], "level="))
1034 1 : nLevel = atoi(papszTokens[i] + 6);
1035 5 : else if (STARTS_WITH_CI(papszTokens[i], "minx="))
1036 : {
1037 1 : bMinXSet = TRUE;
1038 1 : minx = CPLAtof(papszTokens[i] + 5);
1039 : }
1040 4 : else if (STARTS_WITH_CI(papszTokens[i], "miny="))
1041 : {
1042 1 : bMinYSet = TRUE;
1043 1 : miny = CPLAtof(papszTokens[i] + 5);
1044 : }
1045 3 : else if (STARTS_WITH_CI(papszTokens[i], "maxx="))
1046 : {
1047 1 : bMaxXSet = TRUE;
1048 1 : maxx = CPLAtof(papszTokens[i] + 5);
1049 : }
1050 2 : else if (STARTS_WITH_CI(papszTokens[i], "maxy="))
1051 : {
1052 1 : bMaxYSet = TRUE;
1053 1 : maxy = CPLAtof(papszTokens[i] + 5);
1054 : }
1055 1 : else if (STARTS_WITH_CI(papszTokens[i], "bands="))
1056 : {
1057 1 : nReqBands = atoi(papszTokens[i] + 6);
1058 : }
1059 : else
1060 : {
1061 0 : CPLError(CE_Warning, CPLE_AppDefined, "Invalid option : %s",
1062 0 : papszTokens[i]);
1063 : }
1064 : }
1065 17 : CSLDestroy(papszTokens);
1066 : }
1067 :
1068 : /* -------------------------------------------------------------------- */
1069 : /* Open underlying OGR DB */
1070 : /* -------------------------------------------------------------------- */
1071 :
1072 : GDALDatasetH hDS =
1073 54 : RasterliteOpenSQLiteDB(osFileName.c_str(), poOpenInfo->eAccess);
1074 54 : CPLDebug("RASTERLITE", "SQLite DB Open");
1075 :
1076 54 : RasterliteDataset *poDS = nullptr;
1077 :
1078 54 : if (hDS == nullptr)
1079 11 : goto end;
1080 :
1081 43 : if (osTableName.empty())
1082 : {
1083 31 : int nCountSubdataset = 0;
1084 31 : const int nLayers = GDALDatasetGetLayerCount(hDS);
1085 : /* --------------------------------------------------------------------
1086 : */
1087 : /* Add raster layers as subdatasets */
1088 : /* --------------------------------------------------------------------
1089 : */
1090 194 : for (int i = 0; i < nLayers; i++)
1091 : {
1092 163 : OGRLayerH hLyr = GDALDatasetGetLayer(hDS, i);
1093 326 : const std::string osLayerName = OGR_L_GetName(hLyr);
1094 163 : const auto nPosMetadata = osLayerName.find("_metadata");
1095 163 : if (nPosMetadata != std::string::npos)
1096 : {
1097 : const std::string osShortName =
1098 60 : osLayerName.substr(0, nPosMetadata);
1099 :
1100 : const std::string osRasterTableName =
1101 60 : std::string(osShortName).append("_rasters");
1102 :
1103 30 : if (GDALDatasetGetLayerByName(hDS, osRasterTableName.c_str()) !=
1104 : nullptr)
1105 : {
1106 30 : if (poDS == nullptr)
1107 : {
1108 30 : poDS = new RasterliteDataset();
1109 30 : osTableName = osShortName;
1110 : }
1111 :
1112 30 : std::string osSubdatasetName;
1113 30 : if (!STARTS_WITH_CI(poOpenInfo->pszFilename, "RASTERLITE:"))
1114 28 : osSubdatasetName += "RASTERLITE:";
1115 30 : osSubdatasetName += poOpenInfo->pszFilename;
1116 30 : osSubdatasetName += ",table=";
1117 30 : osSubdatasetName += osShortName;
1118 30 : poDS->AddSubDataset(osSubdatasetName.c_str());
1119 :
1120 30 : nCountSubdataset++;
1121 : }
1122 : }
1123 : }
1124 :
1125 31 : if (nCountSubdataset == 0)
1126 : {
1127 1 : goto end;
1128 : }
1129 30 : else if (nCountSubdataset != 1)
1130 : {
1131 0 : poDS->SetDescription(poOpenInfo->pszFilename);
1132 0 : goto end;
1133 : }
1134 :
1135 : /* --------------------------------------------------------------------
1136 : */
1137 : /* If just one subdataset, then open it */
1138 : /* --------------------------------------------------------------------
1139 : */
1140 30 : delete poDS;
1141 30 : poDS = nullptr;
1142 : }
1143 :
1144 : /* -------------------------------------------------------------------- */
1145 : /* Build dataset */
1146 : /* -------------------------------------------------------------------- */
1147 : {
1148 : GDALDataType eDataType;
1149 :
1150 42 : const CPLString osMetadataTableName = osTableName + "_metadata";
1151 :
1152 : OGRLayerH hMetadataLyr =
1153 42 : GDALDatasetGetLayerByName(hDS, osMetadataTableName.c_str());
1154 42 : if (hMetadataLyr == nullptr)
1155 0 : goto end;
1156 :
1157 42 : const CPLString osRasterTableName = osTableName + "_rasters";
1158 :
1159 : OGRLayerH hRasterLyr =
1160 42 : GDALDatasetGetLayerByName(hDS, osRasterTableName.c_str());
1161 42 : if (hRasterLyr == nullptr)
1162 0 : goto end;
1163 :
1164 : /* --------------------------------------------------------------------
1165 : */
1166 : /* Fetch resolutions */
1167 : /* --------------------------------------------------------------------
1168 : */
1169 42 : CPLString osSQL;
1170 : OGRLayerH hSQLLyr;
1171 42 : int nResolutions = 0;
1172 :
1173 : OGRLayerH hRasterPyramidsLyr =
1174 42 : GDALDatasetGetLayerByName(hDS, "raster_pyramids");
1175 42 : if (hRasterPyramidsLyr)
1176 : {
1177 : osSQL.Printf("SELECT pixel_x_size, pixel_y_size "
1178 : "FROM raster_pyramids WHERE table_prefix = '%s' "
1179 : "ORDER BY pixel_x_size ASC",
1180 6 : osTableName.c_str());
1181 :
1182 : hSQLLyr =
1183 6 : GDALDatasetExecuteSQL(hDS, osSQL.c_str(), nullptr, nullptr);
1184 6 : if (hSQLLyr != nullptr)
1185 : {
1186 6 : nResolutions =
1187 6 : static_cast<int>(OGR_L_GetFeatureCount(hSQLLyr, TRUE));
1188 6 : if (nResolutions == 0)
1189 : {
1190 0 : GDALDatasetReleaseResultSet(hDS, hSQLLyr);
1191 0 : hSQLLyr = nullptr;
1192 : }
1193 : }
1194 : }
1195 : else
1196 36 : hSQLLyr = nullptr;
1197 :
1198 42 : if (hSQLLyr == nullptr)
1199 : {
1200 : osSQL.Printf("SELECT DISTINCT(pixel_x_size), pixel_y_size "
1201 : "FROM \"%s_metadata\" WHERE pixel_x_size != 0 "
1202 : "ORDER BY pixel_x_size ASC",
1203 36 : osTableName.c_str());
1204 :
1205 : hSQLLyr =
1206 36 : GDALDatasetExecuteSQL(hDS, osSQL.c_str(), nullptr, nullptr);
1207 36 : if (hSQLLyr == nullptr)
1208 0 : goto end;
1209 :
1210 36 : nResolutions =
1211 36 : static_cast<int>(OGR_L_GetFeatureCount(hSQLLyr, TRUE));
1212 :
1213 36 : if (nResolutions == 0)
1214 : {
1215 0 : GDALDatasetReleaseResultSet(hDS, hSQLLyr);
1216 0 : goto end;
1217 : }
1218 : }
1219 :
1220 : /* --------------------------------------------------------------------
1221 : */
1222 : /* Set dataset attributes */
1223 : /* --------------------------------------------------------------------
1224 : */
1225 :
1226 42 : poDS = new RasterliteDataset();
1227 42 : poDS->SetDescription(poOpenInfo->pszFilename);
1228 42 : poDS->eAccess = poOpenInfo->eAccess;
1229 42 : poDS->osTableName = osTableName;
1230 42 : poDS->osFileName = osFileName;
1231 42 : poDS->hDS = hDS;
1232 :
1233 : /* poDS will release it from now */
1234 42 : hDS = nullptr;
1235 :
1236 : /* --------------------------------------------------------------------
1237 : */
1238 : /* Fetch spatial extent or use the one provided by the user */
1239 : /* --------------------------------------------------------------------
1240 : */
1241 42 : OGREnvelope oEnvelope;
1242 42 : if (bMinXSet && bMinYSet && bMaxXSet && bMaxYSet)
1243 : {
1244 1 : oEnvelope.MinX = minx;
1245 1 : oEnvelope.MinY = miny;
1246 1 : oEnvelope.MaxX = maxx;
1247 1 : oEnvelope.MaxY = maxy;
1248 : }
1249 : else
1250 : {
1251 : CPLConfigOptionSetter oSetter("OGR_SQLITE_EXACT_EXTENT", "YES",
1252 82 : false);
1253 41 : OGR_L_GetExtent(hMetadataLyr, &oEnvelope, TRUE);
1254 : // printf("minx=%.15f miny=%.15f maxx=%.15f maxy=%.15f\n",
1255 : // oEnvelope.MinX, oEnvelope.MinY, oEnvelope.MaxX,
1256 : // oEnvelope.MaxY);
1257 : }
1258 :
1259 : /* --------------------------------------------------------------------
1260 : */
1261 : /* Store resolutions */
1262 : /* --------------------------------------------------------------------
1263 : */
1264 42 : poDS->nResolutions = nResolutions;
1265 42 : poDS->padfXResolutions = reinterpret_cast<double *>(
1266 42 : CPLMalloc(sizeof(double) * poDS->nResolutions));
1267 42 : poDS->padfYResolutions = reinterpret_cast<double *>(
1268 42 : CPLMalloc(sizeof(double) * poDS->nResolutions));
1269 :
1270 : {
1271 : // Add a scope for i.
1272 : OGRFeatureH hFeat;
1273 42 : int i = 0;
1274 93 : while ((hFeat = OGR_L_GetNextFeature(hSQLLyr)) != nullptr)
1275 : {
1276 51 : poDS->padfXResolutions[i] = OGR_F_GetFieldAsDouble(hFeat, 0);
1277 51 : poDS->padfYResolutions[i] = OGR_F_GetFieldAsDouble(hFeat, 1);
1278 :
1279 51 : OGR_F_Destroy(hFeat);
1280 :
1281 : #ifdef RASTERLITE_DEBUG
1282 : printf("[%d] xres=%.15f yres=%.15f\n", i, /*ok*/
1283 : poDS->padfXResolutions[i], poDS->padfYResolutions[i]);
1284 : #endif
1285 :
1286 51 : if (poDS->padfXResolutions[i] <= 0 ||
1287 51 : poDS->padfYResolutions[i] <= 0)
1288 : {
1289 0 : CPLError(CE_Failure, CPLE_NotSupported,
1290 : "res=%d, xres=%.15f, yres=%.15f", i,
1291 0 : poDS->padfXResolutions[i],
1292 0 : poDS->padfYResolutions[i]);
1293 0 : GDALDatasetReleaseResultSet(poDS->hDS, hSQLLyr);
1294 0 : delete poDS;
1295 0 : poDS = nullptr;
1296 0 : goto end;
1297 : }
1298 51 : i++;
1299 : }
1300 : }
1301 :
1302 42 : GDALDatasetReleaseResultSet(poDS->hDS, hSQLLyr);
1303 42 : hSQLLyr = nullptr;
1304 :
1305 : /* --------------------------------------------------------------------
1306 : */
1307 : /* Compute raster size, geotransform and projection */
1308 : /* --------------------------------------------------------------------
1309 : */
1310 42 : const double dfRasterXSize =
1311 42 : (oEnvelope.MaxX - oEnvelope.MinX) / poDS->padfXResolutions[0] + 0.5;
1312 42 : const double dfRasterYSize =
1313 42 : (oEnvelope.MaxY - oEnvelope.MinY) / poDS->padfYResolutions[0] + 0.5;
1314 42 : if (!(dfRasterXSize >= 1 && dfRasterXSize <= INT_MAX) ||
1315 33 : !(dfRasterYSize >= 1 && dfRasterYSize <= INT_MAX))
1316 : {
1317 9 : delete poDS;
1318 9 : poDS = nullptr;
1319 9 : goto end;
1320 : }
1321 33 : poDS->nRasterXSize = static_cast<int>(dfRasterXSize);
1322 33 : poDS->nRasterYSize = static_cast<int>(dfRasterYSize);
1323 :
1324 33 : poDS->bValidGeoTransform = TRUE;
1325 33 : poDS->adfGeoTransform[0] = oEnvelope.MinX;
1326 33 : poDS->adfGeoTransform[1] = poDS->padfXResolutions[0];
1327 33 : poDS->adfGeoTransform[2] = 0;
1328 33 : poDS->adfGeoTransform[3] = oEnvelope.MaxY;
1329 33 : poDS->adfGeoTransform[4] = 0;
1330 33 : poDS->adfGeoTransform[5] = -poDS->padfYResolutions[0];
1331 :
1332 33 : OGRSpatialReferenceH hSRS = OGR_L_GetSpatialRef(hMetadataLyr);
1333 33 : if (hSRS)
1334 : {
1335 16 : poDS->m_oSRS = *(OGRSpatialReference::FromHandle(hSRS));
1336 : }
1337 :
1338 : /* --------------------------------------------------------------------
1339 : */
1340 : /* Get number of bands and block size */
1341 : /* --------------------------------------------------------------------
1342 : */
1343 :
1344 : int nBands;
1345 : int nBlockXSize, nBlockYSize;
1346 33 : if (poDS->GetBlockParams(hRasterLyr, 0, &nBands, &eDataType,
1347 33 : &nBlockXSize, &nBlockYSize) == FALSE)
1348 : {
1349 0 : CPLError(CE_Failure, CPLE_AppDefined,
1350 : "Cannot find block characteristics");
1351 0 : delete poDS;
1352 0 : poDS = nullptr;
1353 0 : goto end;
1354 : }
1355 :
1356 33 : if (eDataType == GDT_Byte && nBands == 1 && nReqBands == 3)
1357 1 : nBands = 3;
1358 32 : else if (nReqBands != 0)
1359 : {
1360 0 : CPLError(CE_Warning, CPLE_NotSupported,
1361 : "Parameters bands=%d ignored", nReqBands);
1362 : }
1363 :
1364 : /* --------------------------------------------------------------------
1365 : */
1366 : /* Add bands */
1367 : /* --------------------------------------------------------------------
1368 : */
1369 :
1370 88 : for (int iBand = 0; iBand < nBands; iBand++)
1371 55 : poDS->SetBand(iBand + 1,
1372 : new RasterliteBand(poDS, iBand + 1, eDataType,
1373 55 : nBlockXSize, nBlockYSize));
1374 :
1375 : /* --------------------------------------------------------------------
1376 : */
1377 : /* Add overview levels as internal datasets */
1378 : /* --------------------------------------------------------------------
1379 : */
1380 33 : if (nResolutions > 1)
1381 : {
1382 6 : poDS->papoOverviews = reinterpret_cast<RasterliteDataset **>(
1383 6 : CPLCalloc(nResolutions - 1, sizeof(RasterliteDataset *)));
1384 15 : for (int nLev = 1; nLev < nResolutions; nLev++)
1385 : {
1386 : int nOvrBands;
1387 : GDALDataType eOvrDataType;
1388 9 : if (poDS->GetBlockParams(hRasterLyr, nLev, &nOvrBands,
1389 : &eOvrDataType, &nBlockXSize,
1390 9 : &nBlockYSize) == FALSE)
1391 : {
1392 0 : CPLError(
1393 : CE_Failure, CPLE_AppDefined,
1394 : "Cannot find block characteristics for overview %d",
1395 : nLev);
1396 0 : delete poDS;
1397 0 : poDS = nullptr;
1398 0 : goto end;
1399 : }
1400 :
1401 9 : if (eDataType == GDT_Byte && nOvrBands == 1 && nReqBands == 3)
1402 0 : nOvrBands = 3;
1403 :
1404 9 : if (nBands != nOvrBands || eDataType != eOvrDataType)
1405 : {
1406 0 : CPLError(CE_Failure, CPLE_AppDefined,
1407 : "Overview %d has not the same number "
1408 : "characteristics as main band",
1409 : nLev);
1410 0 : delete poDS;
1411 0 : poDS = nullptr;
1412 0 : goto end;
1413 : }
1414 :
1415 9 : poDS->papoOverviews[nLev - 1] =
1416 9 : new RasterliteDataset(poDS, nLev);
1417 :
1418 24 : for (int iBand = 0; iBand < nBands; iBand++)
1419 : {
1420 15 : poDS->papoOverviews[nLev - 1]->SetBand(
1421 : iBand + 1, new RasterliteBand(
1422 15 : poDS->papoOverviews[nLev - 1], iBand + 1,
1423 15 : eDataType, nBlockXSize, nBlockYSize));
1424 : }
1425 : }
1426 : }
1427 :
1428 : /* --------------------------------------------------------------------
1429 : */
1430 : /* Select an overview if the user has requested so */
1431 : /* --------------------------------------------------------------------
1432 : */
1433 33 : if (nLevel == 0)
1434 : {
1435 : }
1436 1 : else if (nLevel >= 1 && nLevel <= nResolutions - 1)
1437 : {
1438 1 : poDS->papoOverviews[nLevel - 1]->bMustFree = TRUE;
1439 1 : poDS = poDS->papoOverviews[nLevel - 1];
1440 : }
1441 : else
1442 : {
1443 0 : CPLError(CE_Failure, CPLE_AppDefined,
1444 : "Invalid requested level : %d. Must be >= 0 and <= %d",
1445 : nLevel, nResolutions - 1);
1446 0 : delete poDS;
1447 0 : poDS = nullptr;
1448 : }
1449 : }
1450 :
1451 33 : if (poDS)
1452 : {
1453 : /* --------------------------------------------------------------------
1454 : */
1455 : /* Setup PAM info for this subdatasets */
1456 : /* --------------------------------------------------------------------
1457 : */
1458 33 : poDS->SetPhysicalFilename(osFileName.c_str());
1459 :
1460 66 : CPLString osSubdatasetName;
1461 : osSubdatasetName.Printf("RASTERLITE:%s:table=%s", osFileName.c_str(),
1462 33 : osTableName.c_str());
1463 33 : poDS->SetSubdatasetName(osSubdatasetName.c_str());
1464 33 : poDS->TryLoadXML();
1465 33 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
1466 : }
1467 :
1468 0 : end:
1469 54 : if (hDS)
1470 1 : GDALClose(hDS);
1471 :
1472 54 : return poDS;
1473 : }
1474 :
1475 : /************************************************************************/
1476 : /* GDALRegister_Rasterlite() */
1477 : /************************************************************************/
1478 :
1479 1595 : void GDALRegister_Rasterlite()
1480 :
1481 : {
1482 1595 : if (!GDAL_CHECK_VERSION("Rasterlite driver"))
1483 0 : return;
1484 :
1485 1595 : if (GDALGetDriverByName(DRIVER_NAME) != nullptr)
1486 302 : return;
1487 :
1488 1293 : GDALDriver *poDriver = new GDALDriver();
1489 1293 : RasterliteDriverSetCommonMetadata(poDriver);
1490 :
1491 1293 : poDriver->pfnOpen = RasterliteDataset::Open;
1492 1293 : poDriver->pfnCreateCopy = RasterliteCreateCopy;
1493 1293 : poDriver->pfnDelete = RasterliteDelete;
1494 :
1495 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1496 : }
|
