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.
build overviews on the vrt. As mentioned above, this takes 2+ days to complete for 8 levels. That is unacceptable for my purposes.
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:
Manually create the global pyramids myself. I am wary of recombining them into an .ovr file as I assume that will be tricky.
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.
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.