Displaying huge geotiffs (or vrts) with QGIS?

Question

I am manipulating and processing global rasters at 30m resolution. The total raster sizes are usually [1,440,000 560,000]. I have access to a supercomputer, so I have written code that allows me to break up global rasters into manageable chunks, perform some calculations in parallel, and write them to disk quite rapidly.

I have hit a wall, though, when it comes to displaying results. I usually build a virtual raster of tiles covering the globe and pull it into QGIS. But it's incredibly slow (minutes to load, if it does). And if I try to pan or zoom, it's another many minutes. My first approach to solving this problem was to build overviews using gdaladdo. However, these take forever to build (as in days), which is not conducive to developing algorithms. Here is a list of the things I have tried and why/how they failed.

1. build overviews on the vrt. As mentioned above, this takes 2+ days to complete for 8 levels. That is unacceptable for my purposes.

2. build overviews on the individual tiles, then somehow merge into a vrt that contains the overviews. I am able to build overviews on the tiles rather quickly (supercomputer), but I have not been able to remerge them. I tried:

   2a. gdal_merge on the tiles with overviews, but the overviews were not retained (or at least not recognized by QGIS) in the output tiff.

   2b. gdalbuildvrt on the tiles with overviews, but as above, the overviews were not retained. [This is not correct, see edit.]

   2c. I also tried a hybrid of building overviews for the tiles for levels 1-6 and building levels 7-8 directly on the vrt (basically option 2b) but it is still taking forever for just these two levels. I did some testing and see that the tile overviews are actually used to build to vrt overviews, but it's still on the order of a day to complete the overviews on the vrt.

So I am hoping someone here has some suggestions about where I should go next. Here are some options I'm considering:

1. Manually create the global pyramids myself. I am wary of recombining them into an .ovr file as I assume that will be tricky.

2. Use a mapserver (Geoserver). I know very little about this and am worried that it won't overcome the time hurdles while adding complexity to my process.

3. Split the domain by continents or some other region. I really want to avoid this option.

You might ask "why do you need to view the entire globe at 30m resolution?" One example: I take a mask of water pixels (globally) and skeletonize it to find rivers and perform measurements. My skeletonization algorithm requires a bit of tuning (for branch pruning, removing loops, general cleaning, etc.), and the output is necessarily at 30m. As rivers and landscapes are diverse across the globe, I need to be able to pan around to see the effects of any changes I've implemented.

I have also looked through QGIS to make sure there aren't any settings I could play with to render huge rasters faster, but I didn't see anything. Short of buying SSD drives, I think it's cranking as fast as possible. (My HDDs have I/O of ~250MB/s).

---

I discovered that building overviews on individual tiles, then building a vrt does apparently maintain the overviews--QGIS's "Pyramid" section in the metadata for the file is empty, but in the "Dimensions" section there is an entry for each level of overview (e.g. X 720000, Y 140; X 360000, Y 70, etc.). So I was wrong about 2b.

I also find that if I just pull all the tiles into QGIS, it renders in under a minute, while if I pull in the vrt that references the tiles, it takes >5 minutes (don't know how long exactly since I killed the process).

---

I did some testing on a computer with a SSD, and I found that I could load, display, and render the global vrts (without any overviews) sucessfully and at an acceptable rate. I have ordered a 1TB PCIe SSD in hopes that it will allow me to do the same on my computer. Will update with results.

Answer 1

You seem to have two main concerns: VRT id slow with browsing and it is slow to build global overviews.

---

While I am sure that GDAL VRT used to be slow for me and my MapServer many years ago it may be that the situation has changed. I made a test layer with 10000 aerial images (image/tile size from 10000x10000 to 12000x12000 pixels) and now GDAL VRT is actually faster than native MapServer shapefile index and serves with the test computer 6 tiles (256x256) per second in a simple test with 1 thread where GetMaps hit always the first overview level. Mosaic with 10000 images is still quite a small one and I guess that in my test Linux had the whole VRT file in cache memory. How many images do you have in your VRT?

The following chapter may contain old information, read responsibly:

There is some evidence that VRT is slow when it contains huge number of images. That't because VRT is an index in XML format and it does not support spatial index which leads to full scan of the whole XML file every time. There is nothing you can do for improving that with plain GDAL even there has been some discussion about implementing spatial index for VRT http://osgeo-org.1560.x6.nabble.com/gdal-dev-Don-t-we-have-any-ideas-for-GSoC-2017-td5309810.html.

If you are willing to install new software the easiest workaround could be to use MapServer with tileindex http://www.mapserver.org/optimization/tileindex.html. If you create a tileindex with gdaltindex http://www.gdal.org/gdaltindex.html and create an index for the tileindex as well with shptree http://www.mapserver.org/utilities/shptree.html then MapServer should be able to access very fast all the image files that you have. Create overviews for individual tiles and serve the layer through WMS for QGIS and you have resolved the first part of the problem but not the problem with global overviews. Even if you have created overviews for the individual tiles it will be slow to open thousands of image files for covering a large area and therefore you must limit the number of files by creating overview images which cover larger area. That's what you have already tried to do by building overviews for the hole VRT with gdaladdo.

I do not know any ready made tool in GDAL/MapServer world for creating global pyramids automatically. You could convert tiles from the global VRT into a set of images with bigger pixel size by writing a script that runs gdal_translate http://www.gdal.org/gdal_translate.html with a sliding -prowjin or -srswin. Then you can combine the resulting tiles into a new overview layer with gdalbuildvrt or gdaltindex.

Because you also consider using GeoServer I would recommend to have a loot at gdal_retile script http://www.gdal.org/gdal_retile.html that is written to handle your case. It could be also possible to use the tiles which gdal_retile creates directly as overviews with QGIS by building VRT over them. However, the first problem with slow huge VRT files would remain.

Answer 2

Ok, well I solved both my problems...mainly by buying an NVMe SSD. My disk read/write has gone from 125 MB/s to 1200 MB/s.

Programatically, there are a few things you can do to help your read/write speed. First, consider the blocksize of your tiff. If you are using a striped tiff, when you zoom to a particular region, the GIS software will have to read each complete row of the region, including the portions of the tiff that won't be displayed, in order to display the region. For example, if you zoom into a 256 x 256 pixel region, if you have a striped tiff the software will have to read at least 256 blocks (one per row). If you have a tiled tiff (tiled at 256 x 256), the maximum number of blocks that must be read is 4 (and a minimum of 1). So the first thing you can do is ensure that you are using a tiled tiff (TILED=YES creation option in gdal), and you can set the blocksize to something reasonable (I used 256 x 256 with the gdal creation options BLOCKXSIZE and BLOCKYSIZE.)

Secondly, a hybrid approach to overviews seems to work well. If you can parallelize your operations, you can add overviews to the individual tiles quite quickly, but this will only benefit you for resolutions smaller than your tilesize. I created internal overviews of levels 2 4 8 16 32 and 64 on the individual tiles. Then build a VRT and create overviews of levels 128, 256, and 512 on the VRT (keep in mind that these are for global datasets at 30m resolution--your levels will change depending on number of pixels in your tiff). The total time for creating individual overviews is on the order of minutes (depends on how many threads you can run and how many tiles you have), but creating overviews on the VRT is still on the order of an hour. The runtime improvement over my initial post is due to the SSD and creating fewer levels on the VRT.

Thirdly, you can play with the GDAL_MAX_DATASET_POOL_SIZE option when building vrts as described at the bottom of this page. It sets the maximum number of tiffs to keep in memory at once.

Fourthly, I found that compressing with PACKBITS provides the fastest display times. The files aren't as small as LZW, but it's a tradeoff you might be willing to make.

The result is a VRT that loads rapidly and pans/zooms almost seamlessly.