1

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 with terrain 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:

  1. Create virtual raster
  2. Export new, overlapping tiles (in my case like 100 tiles with about 10% overlap)
  3. Calculate the TPI per tile
  4. Shrink the tiles by the overlap (to get rid of the edge effect)
  5. 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?

5
  • 2
    TIFF (and therefore GeoTIFF) is limited to 4GiB. There exists a BigTIFF spec, but it requires special handling.
    – Vince
    Commented Jan 16 at 14:20
  • Divide and conquer is probably the way to go. But the approach may need to differ depending on the calculation. The TPI is fairly easy to do on a tile basis, as long as you have the surrounding tiles - but still requires some work as you have pointed out. Water flow will be harder, as you need to break the data down by sub catchment... If you get a process for that, I'd be interested in seeing it!! Commented Jan 17 at 5:43
  • You can have look how it's done in lidR package, especially LAScatalog engine: r-lidar.github.io/lidRbook/engine.html. In short - it takes a tile, adds required buffer from surrounding tiles, calculate what's required and returns it to the user. It works for lidars points, but very similar approach might be taken for rasters. Commented Jan 17 at 10:35
  • If you have a hydrography layer with various levels of watersheds, you could use that to create rectangular dems for each of the largest workable watershed areas for flow analysis leaving only the very largest streams that cross watersheds to address, if needed, in a different way.
    – John
    Commented Jan 18 at 15:35
  • Thanks, that would be something to try out (along my idea above). However, it's unlikely I'll find the time to do so. I'm surprised there doesn't seem to be an "out of the box" solution for those problems.
    – Beni
    Commented Jan 19 at 8:52

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.