Input: A large DEM as tiled GeoTIFFs (initially 1'665 tiles, each 500 by 500px, later possibly more).
Goal: To calculate some terrain attributes such as the TPI, waterflows.
Approach for smaller areas: I created a virtual raster and used the TPI from the terrain shading
extension (all in QGIS).
Problem: The algorithm exits after a doing about a quarter, resulting in a single TIFF with a size of about 4 GB (size limit of GeoTIFFs). So the area of interest is too large to be processed at once. When calculating the TPI per tile there are artifacts, where the tiles meet.
Tried:
- Doing the same with R (
saga$ta_morphometry$topographic_position_index_tpi
): It took about 6 hours and completed the job. The exported TIFF is about 2 GB in size. - I used the gdal TPI-algorithm. However, this requires to run
r.filter.stats
prior for good results. It works, with creating TIFF's of around 3 GB. The results are not nearly as good as withterrain shading
. - I used a
whitebox
algorithm for the TPI and they seem to be more data efficient. The resulting tiff was only about 2 GB and the results are good.
However, I'm getting closer to the file-size limit of TIFF's.
Idea:
- Create virtual raster
- Export new, overlapping tiles (in my case like 100 tiles with about 10% overlap)
- Calculate the TPI per tile
- Shrink the tiles by the overlap (to get rid of the edge effect)
- Merge the tiles back to a virtual raster
For flow calculations, I would need to sort the tile(s) with decreasing highest point(s) and start from the highest ones. When calculate adjacent tiles, it would require somehow considering the "outflows" of the previous tile as "inflows".
This would require quite some Python coding. Isn't there a more elegant solution?