7

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.

  • Do you create overviews for VRT file or individual images from VRT file? GDAL allow you to use VRT overview on low resolution (i.e. for 64 x 64, 128 x 128) and use overviews from individual rasters on middle resolution. To create overviews for individual rasters from VRT use gdaladdo on script which will loop throw you files. First of all create individual overviews. Than create overviews for whole VRT. Don't overlap overviews! – Dmitry Baryshnikov Sep 1 '17 at 7:10
  • Feels like the option 2c above and 2+ days for creating the overviews for the whole VRT was considered unacceptable. – user30184 Sep 1 '17 at 7:15
  • Because individual overviews mast be created before. And also it seems to me the overview level was choose wrong. Need to gdaladdo tmp.vrt 2 4 8 16 and gdaladdo individual.tif 2 4 8 16 32 64 ... Topicstarter must calculate the vrt max overview level as: vrt size in pix / image count * 64 pix = max overview for vrt. We created such overviews for north America for Landsat 30 m on desktop PC. The result was created in few days. And rendering of VRT in QGIS was faster than mosaiced geotiff in ArcGIS. – Dmitry Baryshnikov Sep 1 '17 at 7:23
  • 1
    Let's see, they wrote I also tried a hybrid of building overviews for the tiles for levels 1-6 and building levels 7-8 directly on the vrt which feels the same as the order that you recommend and that's naturally the right one. Myself I would not compute 2 4 8 ... overviews for the VRT if individual tiles has them for saving time and disk space. Small ROI would then find overviews from a couple of tiles and that should be fast enough. – user30184 Sep 1 '17 at 7:34
  • This is not true as 2 4 8 of individual file (not tile!) is not the same as 2 4 8 of vrt. For example we have individual image with size 8000 x 8000 and 1000 x 1000 images in vrt. The 64 x 64 pix overview image size will be on 8 level for individual file and for whole vrt it will be on 18 level! And the min level for whole vrt will be 8. So for whole vrt need to use levels from 9 to 18 and for individual files from 2 to 8. – Dmitry Baryshnikov Sep 1 '17 at 7:51
6

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.

  • 1
    I appreciate downvotes but I would like to read also an explanation or a better answer. The MapServer route was my choice when VRT was too slow and it does serve me well with some terabytes of images. – user30184 Sep 1 '17 at 7:25
  • Well, I have been unable to fully complete a global vrt with overviews so I don't know how fast that would render. Without overviews, you are correct that it takes too long to display. I can have as many or as few images in my vrt as I want. I have tried as few as 60 and as many as 4000. Never as many as 10000. I would like to avoid the complexity of a Mapserver-type solution. I did some more research and I think that perhaps the individual tile overviews are indeed retained when I build a vrt--I will test this and post an update on Monday. – Jon Sep 3 '17 at 2:33
4

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.