Skip to main content

Identify flat terrain with above average sun

To identify flat areas in lidar data using GDAL and GRASS in Python, you can use the following steps:

Import the necessary modules​

from osgeo import gdal
import grass.script as gs

Set the GRASS GIS environment​

gisbase = '/usr/local/grass78'
gs.set_gisbase(gisbase)
location = 'location'
mapset = 'mapset'
gs.run_command('g.proj', georef='path/to/georeferenced_file.tif', location=location)

Import the lidar data into GRASS GIS​

gs.run_command('r.in.lidar', input='path/to/lidar_data.las', output='lidar_data', flags='e')

Calculate the slope of the lidar data using the r.slope.aspect module:

gs.run_command('r.slope.aspect', elevation='lidar_data', slope='slope', aspect='aspect')

Identify flat areas by selecting pixels with a slope less than a certain threshold​

gs.mapcalc("flat = if(slope < 0.1, 1, null())")

If you want to also consider sunlight exposure, you can use the r.sun module to calculate solar radiation and sky view factor. For example:

gs.run_command('r.sun', elevation='lidar_data', solar_radiation='solar_radiation', sky_view_factor='sky_view_factor')

You can then use map algebra to select pixels that have both low slope and high solar radiation or sky view factor. For example:

gs.mapcalc("flat_sunny = if(slope < 0.1 && solar_radiation > 500, 1, null())")