1 /**********************************************************************************************
2 Copyright (C) 2014 Oliver Eichler <oliver.eichler@gmx.de>
3
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16
17 **********************************************************************************************/
18
19 #include "CMainWindow.h"
20 #include "helpers/CDraw.h"
21 #include "map/CMapDraw.h"
22 #include "map/CMapVRT.h"
23 #include "units/IUnit.h"
24
25 #include <gdal_priv.h>
26 #include <QtWidgets>
27
28 #define TILELIMIT 2500
29 #define TILESIZEX 64
30 #define TILESIZEY 64
31
32
CMapVRT(const QString & filename,CMapDraw * parent)33 CMapVRT::CMapVRT(const QString& filename, CMapDraw* parent)
34 : IMap(eFeatVisibility, parent)
35 , filename(filename)
36 {
37 qDebug() << "------------------------------";
38 qDebug() << "VRT: try to open" << filename;
39
40 dataset = (GDALDataset*)GDALOpen(filename.toUtf8(), GA_ReadOnly);
41
42 if(nullptr == dataset)
43 {
44 QMessageBox::warning(CMainWindow::getBestWidgetForParent(), tr("Error..."), tr("Failed to load file: %1").arg(filename));
45 return;
46 }
47
48 // ------- setup color table ---------
49 rasterBandCount = dataset->GetRasterCount();
50 if(rasterBandCount == 1)
51 {
52 GDALRasterBand* pBand = dataset->GetRasterBand(1);
53
54 if(nullptr == pBand)
55 {
56 GDALClose(dataset);
57 dataset = nullptr;
58 QMessageBox::warning(CMainWindow::getBestWidgetForParent(), tr("Error..."), tr("Failed to load file: %1").arg(filename));
59 return;
60 }
61
62 if(pBand->GetColorInterpretation() == GCI_PaletteIndex )
63 {
64 GDALColorTable* pct = pBand->GetColorTable();
65 for(int i = 0; i < pct->GetColorEntryCount(); ++i)
66 {
67 const GDALColorEntry& e = *pct->GetColorEntry(i);
68 colortable << qRgba(e.c1, e.c2, e.c3, e.c4);
69 }
70 }
71 else if(pBand->GetColorInterpretation() == GCI_GrayIndex )
72 {
73 for(int i = 0; i < 256; ++i)
74 {
75 colortable << qRgba(i, i, i, 255);
76 }
77 }
78 else
79 {
80 GDALClose(dataset);
81 dataset = nullptr;
82 QMessageBox::warning(CMainWindow::getBestWidgetForParent(), tr("Error..."), tr("File must be 8 bit palette or gray indexed."));
83 return;
84 }
85
86 int success = 0;
87 qreal idx = pBand->GetNoDataValue(&success);
88 if(success)
89 {
90 if((idx >= 0) && (idx < colortable.size()))
91 {
92 QColor tmp(colortable[idx]);
93 tmp.setAlpha(0);
94 colortable[idx] = tmp.rgba();
95 }
96 else
97 {
98 qDebug() << "Index for no data value is out of bound";
99 return;
100 }
101 }
102 }
103
104 if(dataset->GetRasterCount() > 0)
105 {
106 hasOverviews = dataset->GetRasterBand(1)->GetOverviewCount() != 0;
107 }
108
109 // if the master VRT does not return a positive overview feedback
110 // test all files combined by the VRT to have overviews.
111 if(!hasOverviews)
112 {
113 qDebug() << "extended test for overviews";
114 hasOverviews = testForOverviews(filename);
115 }
116 qDebug() << "has overviews" << hasOverviews;
117
118
119 // ------- setup projection ---------------
120 proj.init(dataset->GetProjectionRef(), "EPSG:4326");
121
122 if(!proj.isValid())
123 {
124 delete dataset;
125 dataset = nullptr;
126 QMessageBox::warning(CMainWindow::getBestWidgetForParent(), tr("Error..."), tr("No georeference information found."));
127 return;
128 }
129
130 xsize_px = dataset->GetRasterXSize();
131 ysize_px = dataset->GetRasterYSize();
132
133
134 qreal adfGeoTransform[6];
135 dataset->GetGeoTransform( adfGeoTransform );
136
137 xscale = adfGeoTransform[1];
138 yscale = adfGeoTransform[5];
139 xrot = adfGeoTransform[4];
140 yrot = adfGeoTransform[2];
141
142 trFwd.translate(adfGeoTransform[0], adfGeoTransform[3]);
143 trFwd.scale(adfGeoTransform[1], adfGeoTransform[5]);
144
145 if(adfGeoTransform[4] != 0.0)
146 {
147 trFwd.rotate(qAtan(adfGeoTransform[2] / adfGeoTransform[4]));
148 }
149
150 if(proj.isSrcLatLong())
151 {
152 // convert to RAD to match internal notations
153 trFwd = trFwd * DEG_TO_RAD;
154 }
155
156 trInv = trFwd.inverted();
157
158 ref1 = trFwd.map(QPointF(0, 0));
159 ref2 = trFwd.map(QPointF(xsize_px, 0));
160 ref3 = trFwd.map(QPointF(xsize_px, ysize_px));
161 ref4 = trFwd.map(QPointF(0, ysize_px));
162
163 qDebug() << "FF" << trFwd;
164 qDebug() << "RR" << trInv;
165
166 isActivated = true;
167 }
168
~CMapVRT()169 CMapVRT::~CMapVRT()
170 {
171 GDALClose(dataset);
172 }
173
testForOverviews(const QString & filename)174 bool CMapVRT::testForOverviews(const QString& filename)
175 {
176 QFile file(filename);
177 if (!file.open(QIODevice::ReadOnly))
178 {
179 qDebug() << "Failed to open" << filename;
180 return false;
181 }
182
183 QDomDocument xml;
184 QString msg;
185 int line;
186 int column;
187 if (!xml.setContent(&file, false, &msg, &line, &column))
188 {
189 file.close();
190 throw tr("Failed to read: %1\nline %2, column %3:\n %4").arg(filename).arg(line).arg(column).arg(msg);
191 qDebug() << "Failed to read:" << filename << endl
192 << "line" << line << ", column" << column << endl
193 << msg;
194 return false;
195 }
196 file.close();
197
198 QSet<QString> files;
199 QDir basePath(QFileInfo(filename).absoluteDir());
200 const QDomElement& xmlVrt = xml.documentElement();
201
202 {
203 const QDomNodeList& xmlComplexSources = xmlVrt.elementsByTagName("ComplexSource");
204 const int N = xmlComplexSources.count();
205 for(int n = 0; n < N; ++n)
206 {
207 const QDomNode& xmlComplexSource = xmlComplexSources.item(n);
208 const QDomNode& xmlSourceFilename = xmlComplexSource.namedItem("SourceFilename");
209 const QDomNamedNodeMap& attr = xmlSourceFilename.attributes();
210 QString subFilename = xmlSourceFilename.toElement().text();
211
212 if(attr.contains("relativeToVRT") && (attr.namedItem("relativeToVRT").nodeValue() == "1"))
213 {
214 subFilename = basePath.absoluteFilePath(subFilename);
215 }
216
217 files << subFilename;
218 }
219 }
220
221 {
222 const QDomNodeList& xmlSimpleSources = xmlVrt.elementsByTagName("SimpleSource");
223 const int N = xmlSimpleSources.count();
224 for(int n = 0; n < N; ++n)
225 {
226 const QDomNode& xmlSimpleSource = xmlSimpleSources.item(n);
227 const QDomNode& xmlSourceFilename = xmlSimpleSource.namedItem("SourceFilename");
228 const QDomNamedNodeMap& attr = xmlSourceFilename.attributes();
229 QString subFilename = xmlSourceFilename.toElement().text();
230
231 if(attr.contains("relativeToVRT") && (attr.namedItem("relativeToVRT").nodeValue() == "1"))
232 {
233 subFilename = basePath.absoluteFilePath(subFilename);
234 }
235
236 files << subFilename;
237 }
238 }
239
240 if(files.isEmpty())
241 {
242 return false;
243 }
244
245 for(const QString& file : qAsConst(files))
246 {
247 using pGDALDataset = QSharedPointer<GDALDataset>;
248 pGDALDataset _dataset = pGDALDataset((GDALDataset*)GDALOpen(file.toUtf8(), GA_ReadOnly), GDALClose);
249 // _dataset will be destroyed automatically by shared pointer.
250
251 if(_dataset == nullptr)
252 {
253 return false;
254 }
255
256 if(_dataset->GetRasterCount() > 0)
257 {
258 if(_dataset->GetRasterBand(1)->GetOverviewCount() == 0)
259 {
260 return false;
261 }
262 }
263 else
264 {
265 return false;
266 }
267 }
268
269
270 return true;
271 }
272
draw(IDrawContext::buffer_t & buf)273 void CMapVRT::draw(IDrawContext::buffer_t& buf) /* override */
274 {
275 if(map->needsRedraw())
276 {
277 return;
278 }
279
280 QPointF bufferScale = buf.scale * buf.zoomFactor;
281
282 // calculate bounding box;
283 QPointF pt1 = ref1;
284 QPointF pt2 = ref2;
285 QPointF pt3 = ref3;
286 QPointF pt4 = ref4;
287
288 proj.transform(pt1, PJ_FWD);
289 proj.transform(pt2, PJ_FWD);
290 proj.transform(pt3, PJ_FWD);
291 proj.transform(pt4, PJ_FWD);
292
293 QPolygonF boundingBox;
294 boundingBox << pt1 << pt2 << pt3 << pt4;
295 map->convertRad2Px(boundingBox);
296
297 // get pixel offset of top left buffer corner
298 QPointF pp = buf.ref1;
299 map->convertRad2Px(pp);
300
301 // calculate area to read from file
302 pt1 = buf.ref1;
303 pt2 = buf.ref2;
304 pt3 = buf.ref3;
305 pt4 = buf.ref4;
306
307 proj.transform(pt1, PJ_INV);
308 proj.transform(pt2, PJ_INV);
309 proj.transform(pt3, PJ_INV);
310 proj.transform(pt4, PJ_INV);
311
312 pt1 = trInv.map(pt1);
313 pt2 = trInv.map(pt2);
314 pt3 = trInv.map(pt3);
315 pt4 = trInv.map(pt4);
316
317 qreal left, right, top, bottom;
318 left = pt1.x() < pt4.x() ? pt1.x() : pt4.x();
319 right = pt2.x() > pt3.x() ? pt2.x() : pt3.x();
320 top = pt1.y() < pt2.y() ? pt1.y() : pt2.y();
321 bottom = pt4.y() > pt3.y() ? pt4.y() : pt3.y();
322
323 if(left < 0)
324 {
325 left = 0;
326 }
327 if(left > xsize_px)
328 {
329 left = xsize_px;
330 }
331
332 if(top < 0)
333 {
334 top = 0;
335 }
336 if(top > ysize_px)
337 {
338 top = ysize_px;
339 }
340
341 if(right > xsize_px)
342 {
343 right = xsize_px;
344 }
345 if(right < 0)
346 {
347 right = 0;
348 }
349
350 if(bottom > ysize_px)
351 {
352 bottom = ysize_px;
353 }
354 if(bottom < 0)
355 {
356 bottom = 0;
357 }
358
359 qint32 imgw = TILESIZEX;
360 qint32 imgh = TILESIZEY;
361 qint32 dx = imgw;
362 qint32 dy = imgh;
363
364
365 // estimate number of tiles and use it as a limit if no
366 // user defined limit is given
367 qreal nTiles = ((right - left) * (bottom - top) / (dx * dy));
368 if(hasOverviews)
369 {
370 // if there are overviews tiles can be reduced by reading
371 // with a scale factor from file. Increase amount of pixel
372 // read until tile limit is met.
373 while(nTiles > TILELIMIT)
374 {
375 dx *= 2;
376 dy *= 2;
377 nTiles /= 4;
378 }
379 }
380 else
381 {
382 nTiles = getMaxScale() == NOFLOAT ? nTiles : 0;
383 }
384
385 // start to draw the map
386 QPainter p(&buf.image);
387 USE_ANTI_ALIASING(p, true);
388 p.setOpacity(getOpacity() / 100.0);
389 p.translate(-pp);
390
391
392 // qDebug() << imgw << dx << nTiles;
393 // limit number of tiles to keep performance
394 if(!isOutOfScale(bufferScale) && (nTiles < TILELIMIT))
395 {
396 for(qint32 y = top; y < bottom; y += dy)
397 {
398 if(map->needsRedraw())
399 {
400 break;
401 }
402
403 for(qint32 x = left; x < right; x += dx)
404 {
405 if(map->needsRedraw())
406 {
407 break;
408 }
409
410 // read tile from file
411 CPLErr err = CE_Failure;
412
413 // reduce tile size at the border of the file
414 qreal dx_used = dx;
415 qreal dy_used = dy;
416 qreal imgw_used = imgw;
417 qreal imgh_used = imgh;
418
419 if((x + dx) > xsize_px)
420 {
421 dx_used = xsize_px - x;
422 imgw_used = qRound(imgw * dx_used / dx) & 0xFFFFFFFC;
423 }
424 if((y + dy) > ysize_px)
425 {
426 dy_used = ysize_px - y;
427 imgh_used = imgh * dy_used / dy;
428 }
429
430 dx_used = qFloor(dx_used);
431 dy_used = qFloor(dy_used);
432 imgw_used = qRound(imgw_used);
433 imgh_used = qRound(imgh_used);
434
435 if(imgw_used < 1 || imgh_used < 1)
436 {
437 continue;
438 }
439
440 QImage img;
441 if(rasterBandCount == 1)
442 {
443 GDALRasterBand* pBand;
444 pBand = dataset->GetRasterBand(1);
445
446 img = QImage(QSize(imgw_used, imgh_used), QImage::Format_Indexed8);
447 img.setColorTable(colortable);
448
449 err = pBand->RasterIO(GF_Read
450 , x, y
451 , dx_used, dy_used
452 , img.bits()
453 , imgw_used, imgh_used
454 , GDT_Byte, 0, 0);
455 }
456 else
457 {
458 img = QImage(imgw_used, imgh_used, QImage::Format_ARGB32);
459 img.fill(qRgba(255, 255, 255, 255));
460
461 QVector<quint8> buffer(imgw_used* imgh_used);
462
463 QRgb testPix = qRgba(GCI_RedBand, GCI_GreenBand, GCI_BlueBand, GCI_AlphaBand);
464
465 for(int b = 1; b <= rasterBandCount; ++b)
466 {
467 GDALRasterBand* pBand;
468 pBand = dataset->GetRasterBand(b);
469
470 err = pBand->RasterIO(GF_Read
471 , x, y
472 , dx_used, dy_used
473 , buffer.data()
474 , imgw_used, imgh_used
475 , GDT_Byte, 0, 0);
476
477 if(!err)
478 {
479 int pbandColour = pBand->GetColorInterpretation();
480 unsigned int offset;
481
482 for (offset = 0; offset < sizeof(testPix) && *(((quint8*)&testPix) + offset) != pbandColour; offset++)
483 {
484 }
485 if(offset < sizeof(testPix))
486 {
487 quint8* pTar = img.bits() + offset;
488 quint8* pSrc = buffer.data();
489 const int size = buffer.size();
490
491 for(int i = 0; i < size; ++i)
492 {
493 *pTar = *pSrc;
494 pTar += sizeof(testPix);
495 pSrc += 1;
496 }
497 }
498 }
499 }
500 }
501
502 if(err)
503 {
504 continue;
505 }
506
507
508 QPolygonF l;
509 l << QPointF(x, y) << QPointF(x + dx_used, y) << QPointF(x + dx_used, y + dy_used) << QPointF(x, y + dy_used);
510 l = trFwd.map(l);
511
512 proj.transform(l, PJ_FWD);
513
514 drawTile(img, l, p);
515 }
516 }
517 }
518
519 p.setPen(Qt::black);
520 p.setBrush(Qt::NoBrush);
521 p.drawPolygon(boundingBox);
522 }
523
524