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