[gdal-dev] Processing huge DEM dataset

[Even Rouault]

Le jeudi 23 avril 2015 14:22:30, Paul Ramsey a écrit :
> I have in the past, with other tool sets, not GDAL, approached this by
> building out padded tiles as the first step. So for each tile, merge it
> with it’s neighbors, then clip out the middle so you get a somewhat larger
> tile. Give it a nice thick buffer.
> 
> Now all your tiles overlap. Process them all individually. Thanks to the
> generous overlap they should, on the boundaries, come up w/ the same
> answers. Once you have final results (hillshades, contours, etc) clip them
> to the real tile boundary. In the case of contours you may need to finally
> snap the ends together, but it should not be a huge snap, just a tiny one.
> 

+1 for Paul's above workflow. Some answers to your questions below

> --
> http://postgis.net
> http://cleverelephant.ca
> 
> On April 23, 2015 at 3:01:45 AM, Marcos Dione (mdione at grulic.org.ar) wrote:
> > I have SRTM's DEM 1x1 degree 30m resolution tiles for the whole Europe
> > and I'm trying to generate several raster images based on
> > that (elevation coloring, slopeshade and hillshade), but I'm not sure
> > about the right approach to do it for that amount of data.
> > 
> > The simplest approach is to stitch the DEMs and then process, but that
> > takes ages, specially if I try to use uncompressed, tiled
> > GeoTIFFs as output. The stitching can even be done using a virtual file,
> > which saves space.
> > 
> > If I process each tile individually, and then build a virtual file on
> > top, I get shades on the edges of tiles. This shade is due
> > to the tile ending and the shading algorithm assuming there's a 0
> > elevation point right to it. So, question A) is that so?

Most gdaldem algorithms (except color-relief) need to compute some form of 
gradient (a 3x3 window around the pixel being computed), so you have edge 
effects. By default, they put a nodata value on the edges.
If you specify -compute_edges, then they will interpolate extra values from 
the ones available so that the edge pixels can be computed. You could still 
see some discontinuity if the prediction isn't that great.

> > 
> > I think that getting the shade in the output is due to the algorithm for
> > finding a pixel uses the first tile that has it.
> > Question B) is that so?
> > 
> > If so, C) could I simply avoid this by generating another vrt file that
> > lists each tile as having a bbox of only the 1x1 degree
> > instead of the 1x1 degree plus an extra pixel border? If I get the time,
> > I'll try this this afternoon (I just thought of it).

For each tile, you could have a VRT of 3x3 tiles. Let's say that a tile if NxN 
pixel, then the VRT would be (N+2)*(N+2). And you would extract NxN pixels 
from the output of gdaldem on that VRT, keeping NxN pixels only. Well, this is 
basically Paul's approach using VRT to do the buffer.

> > 
> > All this I can do more or less with the gdal command line tools, without
> > much programming. Then comes a more programmatic way:
> > either use gdaldem or use the GDAL API to process each tile, then cut the
> > 1x1 degree image and save that; then stitch them/build a
> > vrt file on top. As you can see, this is what I've been avoiding to do :)
> > 
> > Finally, I would also like to generate contour lines for this. So far I
> > managed to generate them for 5x5 tiles with 90m
> > resolution, then I import them in postgis. When I render them, on the
> > edges of such tiles I see the lines from one crossing the
> > others, looking ugly. For instance:
> > 
> > http://grulicueva.homenet.org/~mdione/Elevation/#14/45.0000/15.0000
> > 
> > I tried used a stitched file for the whole region but I ran out of memory
> > with gdal_contour. Again, this was with 90m resolution
> > tiles; now I have 30m, which means 9 times as much data. D) How could I
> > properly process all that?
> > 
> > Thanks in advance for any ideas,
> > 
> > -- Marcos.
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com