Calculate focal range on XYZ DEM in QGIS or Spatialite

Question

I have a XYZ DEM of a fairly mountainous terrain and would like to identify flat areas using QGIS. How can calculate focal range in QGIS? Or would there be a way to do this in Spatialite?

Ideally, the output would be a shapefile of polygons of flat areas or a set of points representing centers of flat circles of certain radii.

Answer 1

This is doable with QGIS and SAGA (either standalone SAGA, or via the Processing panel in QGIS)

I assume by 'XZY DEM' you mean an ascii grid.

If the points are not regularly gridded (e.g. from UAS survey or LIDAR), you can also do this by importing into SAGA as a Point Cloud from an .xyz text file, then converting the point cloud into raster.

If it's an ascii gridded XYZ file, or a normal raster, you can import it using File > GDAL/OGR > Import Raster.

You can find flat areas using Terrain > Preprocessing > Flat Detection.

That will create two rasters;

• one for flat areas only,
• one for non-flat areas.

Depending on what you want, you can choose to set the pixel values for the flat raster to be set to either the elevation of each flat area ('elevation') or a unique ID ('enumeration'). Everything else is set to NODATA.

The screenshot below shows the unique ID for lochs in central Scotland. As far as I can tell, flat areas are numbered left-to-right, bottom-to-top, hence the colour gradient.

You can then convert these to polygons, either in SAGA or in QGIS, and use Simplify Geometries if necessary to remove the pixel artifacts along the edges.

EDIT:

To get an XYZ file into SAGA as a point cloud, use File > Shapes > Import > Point Cloud from Text File. I find the point clouds can be tricky to work with if RAM is tight, so I sometimes use a python script to randomly sample (say) 1% or 10% of points

Next, Shapes > Point Clouds > Conversion > Point Cloud to Grid You probably want the defaults (only Z, first value) but check to see if one of the other options is more apt for your needs (if more than one point maps to a pixel, you can choose the maximum, minimum, average etc).

You can also choose the pixel size. The files I work with have coords in UTM30N (meters), so the cell size is the number of meters per pixel.

Finally, you'll have a raster. Well, actually two rasters..

• MyGrid [Z]
• MyGrid [Points per cell]

The first one has your DEM, the second shows how many points fell into each pixel.

You may well have missing data - pixels where no samples were found (especially if you only work with a subset of the data).

To fix this, you can use either

• Grid > Gaps > Close One cell gaps (fast, but only works with single pixel voids)
• Grid > Gaps > Close Gaps (slow, but works with larger voids)