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