Viewshed analysis in GRASS and Jupyter ¶
This is a project where I calculate viewsheds using GRASS and Jupyter. This notebook is based on the excellent FOSS4G 2022 workshop.
Prepare environment¶
Load dependencies:
import os
import sys
import json
import subprocess
import psutil
import pandas as pd
import matplotlib.pyplot as plt
import time
from pathlib import Path
from tqdm import tqdm
from typing import List, Tuple, Dict, Optional
from IPython.display import clear_output, display, HTML
from multiprocessing import Pool, cpu_count
Are we running in 64bit?
import sys
if sys.maxsize > 2**32:
print('64-bit')
else:
print('32-bit')
How much memory? (MB)
TOTAL_MEMORY_MB = int(psutil.virtual_memory().total / 1024 / 1024)
TOTAL_MEMORY_MB
How many CPUs?
PROC = cpu_count()
PROC
- Ask GRASS GIS where its Python packages are.
- Import the GRASS GIS packages we need.
- Start GRASS Session
sys.path.append(
subprocess.check_output(
["grass", "--config", "python_path"], text=True).strip()
)
import grass.script as gs
import grass.jupyter as gj
from grass.exceptions import CalledModuleError
Activate autoreload of changed python files:
%load_ext autoreload
%autoreload 2
Parameters¶
Define initial parameters that affect processing
MAX_DISTANCE = 10000 # Max distance of viewshed calculation (in Meters)
HA_FILTER_NOISE = 5 # minimal area unit for viewshed noise removal
Current working Path
Path.cwd()
DATA = Path.cwd().parents[0] / "grassdata"
DATA.mkdir(exist_ok=True)
(Path.cwd().parents[0] / "notebooks").mkdir(exist_ok=True)
(Path.cwd().parents[0] / "py").mkdir(exist_ok=True)
Notebook environment¶
By default all cells are running Python:
import sys
v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")
Cells can also run bash using IPython magic. Later, we will use this to show how GRASS can be used from shell.
%%bash
pwd
or we can use ! to run individual lines in bash:
!pwd
Create grass project (a Location)¶
In this first part, we will demonstrate starting GRASS GIS, creating new project and basic data visualization.
Set WORK_FOLDER in python, export to bash/env variable.
PROJECT_NAME = "viewshed"
WORK_FOLDER = "/home/jovyan/work/viewshed_analysis/grassdata/"
PROJECT_FOLDER = f'{WORK_FOLDER}{PROJECT_NAME}/'
os.environ['PROJECT_FOLDER'] = PROJECT_FOLDER
First, we create new empty location (project) called viewshed.
Flag c stands for creating new location and e will exit the command after creating the location. See manual for more examples.
Check if folder exists, otherwise create location.
We use 25833 EPSG projection, which is suitable for Germany and matches out input digital elevation model projection.
%%bash
[[ -d $PROJECT_FOLDER ]] || grass -c EPSG:25833 -e $PROJECT_FOLDER
Start GRASS GIS session in the newly created location.
session = gj.init(WORK_FOLDER, PROJECT_NAME, "PERMANENT")
Change working directory into project folder
os.chdir(f'{WORK_FOLDER}{PROJECT_NAME}')
Import data¶
We will import prepared digitial surface model (DSM). This data was exported from ArcPro ESRI to a Geo-TIFF. The data CRS matches the CRS of the viewshed location, so we don't need to reproject it.
%%time
%%capture output
%%bash
r.import input=dsm.tif output=dsm resample=bilinear
print(output.stderr)
List available data in our project:
- this may list results from previous runs.
%%bash
g.list type=raster,vector -m -t
Data visualization¶
We will call d.rast/d.vect modules using the new GRASS Jupyter API. The Map class creates and displays GRASS maps as static PNG images.
First let's display bare ground elevation:
# Create Map instance
elevation_map = gj.Map()
# Add a dsm raster to the map
elevation_map.d_rast(map="dsm")
# Display map
elevation_map.show()
Load viewpoints¶
We will first derive viewpoints that will be used for the viewshed calculation.
- Load points (Pot Wea 0 Test)
- Create new vector of points dataset
%%capture output
%%bash
v.import input="POT_WEA/pot_WEA_0.shp" output=pot_wea_0
print(output.stderr)
Visualize the points with InteractiveMap with OSM tiles (see other tile options):
gs.run_command("g.region", raster="dsm", save="zoom_region", grow=250)
dsm_map = gj.InteractiveMap(tiles="OpenStreetMap", saved_region="zoom_region")
dsm_map.add_vector("pot_wea_0")
dsm_map.show()
Next part of analysis is raster-based, so we need to make sure we set computational region as we need. Specifically, we set it to match the DSM:
gs.run_command("g.region", raster="dsm")
Viewshed computation¶
To get the cumulative viewshed, we will compute viewsheds from all the viewpoints.
First, we get the list coordinates of the viewpoints and height data (this is a wind turbines dataset):
df = pd.DataFrame(json.loads(gs.read_command("v.db.select", map="pot_wea_0", columns="cat,Gesamthoeh", layer=1, format="json"))["records"])
df
viewpoints = gs.read_command(
'v.out.ascii', input='pot_wea_0',
separator='comma', layer=1).strip().splitlines()
viewpoints = [p.split(",") for p in viewpoints]
[print(p) for p in viewpoints]
viewpoints
We will now compute the viewshed from each viewpoint in a loop. We set max distance of 10km (MAX_DISTANCE parameter). Each viewshed will be named viewshed_cat. We consider the curvature of earth. We assign max 25% of memory to the computation, which is about 32GB on this VM (128GB total).
Configure r.viewshed Paremeters.
vkwargs = {
"input":"dsm", # digital surface model to be used
"max_distance":MAX_DISTANCE, # 10km distance (max) viewshed area radius
"flags":"bc", # compute visible area, consider curvature of earth
"memory":TOTAL_MEMORY_MB/4, # allocate max 25% of max memory per computation
"target_elevation": 1.8, # Offset for target elevation above the ground
}
%%time
maps = []
for x, y, height, cat in tqdm(viewpoints):
name = f"viewshed_{cat}"
gs.run_command(
"r.viewshed", output=name, observer_elevation=height,
coordinates=(x, y), **vkwargs)
maps.append(name)
Since these are independent runs, we can easily parallelize the r.viewshed calls using Python multiprocessing.
We define a function that computes the viewshed and returns the name of the output or None in case of error:
def viewshed(point) -> str:
"""Calculate viewshed for point and return (name)"""
x, y, height, cat = point
# ghoeh = df["Gesamthoeh"][int(cat)-1]
x, y = float(x), float(y)
name = f"viewshed_{cat}"
try:
gs.run_command(
"r.viewshed", output=name, observer_elevation=height,
coordinates=(x, y), **vkwargs)
return f"viewshed_{cat}"
except CalledModuleError:
return None
Run with the total number of CPUs (threads) available (PROC)
%%time
with Pool(processes=PROC) as pool:
maps = pool.map_async(viewshed, viewpoints).get()
print(maps)
One trick to speedup viewshed computation is to limit the computation only to the actual area given by the maxdistance option. To do that we will locally modify the computational region and pass the environment to the module directly. The current computational region won't be affected.
def viewshed(point, max_distance: int = None):
"""Calculate viewshed for point, limit to region, store raster, and return the name"""
x, y, height, cat = point
# ghoeh = df["Gesamthoeh"][int(cat)-1]
x, y = float(x), float(y)
if max_distance is None:
max_distance = MAX_DISTANCE
# copy current environment
env = os.environ.copy()
# set GRASS_REGION variable using region_env function
env["GRASS_REGION"] = gs.region_env(
align="dsm",
e=x + max_distance,
w=x - max_distance,
n=y + max_distance,
s=y - max_distance)
grp_str = ""
name = f"viewshed_{cat}"
try:
gs.run_command(
"r.viewshed", output=name, observer_elevation=height,
coordinates=(x, y), env=env, **vkwargs)
return f"viewshed_{cat}"
except CalledModuleError:
return None
Run with the number of CPUs available
%%time
with Pool(processes=PROC) as pool:
maps = pool.map_async(viewshed, viewpoints).get()
print(maps)
print(f"Viewshed num cells: {gs.raster_info(maps[0])['cells']}")
print(f"DSM num cells: {gs.raster_info('dsm')['cells']}")
preview viewshed 1 of 9:
gs.run_command("g.region", raster=maps, save="zoom_region", grow=250)
tileset = "Esri.WorldImagery" # "OpenStreetMap"
cumulative_map = gj.InteractiveMap(tiles=tileset, saved_region="zoom_region")
cumulative_map.add_raster("viewshed_1", opacity=0.8)
cumulative_map.add_layer_control(position="bottomright")
cumulative_map.add_vector("pot_wea_0")
cumulative_map.show()
Cumulative viewshed¶
We can compute the cumulative viewshed, which aggregates viewsheds from multiple viewpoints. In this way you can e.g., identify the most frequently visible areas for the wind turbines dataset.
Since our viewshed rasters are binary (0 invisible, 1 visible), we will use r.series method sum. As a result, we get a cumulative viewshed raster with the number of overlays encoded. Then we replace zeros with no data using r.null and set a new color ramp:
gs.run_command("r.series", input=maps, output="cumulative_viewshed", method="sum")
gs.run_command("r.null", map="cumulative_viewshed", setnull=0)
gs.run_command("r.colors", map="cumulative_viewshed", color="plasma")
Let's visualize the results:
- Interactive Folium map
- Select the Tile provider
- zoom in to data using
saved_region, leave a250border
gs.run_command("g.region", raster=maps, save="zoom_region", grow=250)
tileset = "Esri.WorldImagery" # "OpenStreetMap"
cumulative_map = gj.InteractiveMap(tiles=tileset, saved_region="zoom_region")
cumulative_map.add_raster("cumulative_viewshed", opacity=0.8)
cumulative_map.add_layer_control(position="bottomright")
cumulative_map.add_vector("pot_wea_0")
cumulative_map.show()
Set non-null to 1 (binary)
# gs.run_command("g.region", save="zoom_region")
gs.mapcalc('cumulative_viewshed_bin = if(isnull(cumulative_viewshed), cumulative_viewshed, 1)')
Cleanup Map based on 5 ha minimum area:
gs.run_command(
"r.reclass.area", input="cumulative_viewshed_bin", output="cumulative_viewshed_cleaned", mode="greater",
value=HA_FILTER_NOISE)
Visualize as static png with the dsm in background.
viewshed_map = gj.Map(saved_region="zoom_region")
viewshed_map.d_rast(map="dsm")
viewshed_map.d_rast(map="cumulative_viewshed_cleaned")
viewshed_map.show()
Remove stop areas¶
Import stop areas
%%time
%%capture output
%%bash
v.import input="stop_areas/sichtverstellende_Elemente.shp" output=stop_areas
Convert to raster
%%time
%%capture output
%%bash
v.to.rast input=stop_areas output=stop_areas_raster use=val val=1
print(output.stderr)
gs.mapcalc(
'cumulative_viewshed_stopareas = if(isnull(stop_areas_raster), cumulative_viewshed_cleaned, 0)')
display_raster_map(raster_layer="cumulative_viewshed_stopareas")
gs.run_command("r.null", map="cumulative_viewshed_stopareas", setnull=0)
display_raster_map(raster_layer="cumulative_viewshed_stopareas")
Export result¶
Export map
viewshed_map.save("viewshed0.png")
Write to GeoTiff
gs.run_command("g.region", region="zoom_region")
gs.run_command(
"r.out.gdal",
flags="fmt",
input="cumulative_viewshed_cleaned",
output=f"{PROJECT_FOLDER}viewshed0.tif",
format="GTiff",
overviews="4")
Cumulative viewshed with r.viewshed.cva¶
ToDo:
- See r.viewshed.cva
- it seems like different observation heights are not possible
Loop through all points¶
Loop through complete wind turbines point dataset
Let's create a few methods first and then we will loop through all *.shp and create individual viewsheds.
def get_viewpoints(shp_name: str):
"""Get group of viewpoints from a given shapefile"""
layer_name = shp_name[:-4]
gs.run_command("v.import", input=f"POT_WEA/{shp_name}", output="pot_wea")
# v.import input=f"POT_WEA/{shp_name}" output=layer_name
viewpoints = gs.read_command(
'v.out.ascii', input="pot_wea",
separator='comma', layer=1).strip().splitlines()
viewpoints = [p.split(",") for p in viewpoints]
return viewpoints
def create_viewshed_group(viewpoints):
"""Create a viewshed for a given group of points"""
with Pool(processes=PROC) as pool:
maps = pool.map_async(viewshed, viewpoints).get()
return maps
def cumulative_viewshed(maps, output_name: str):
"""Calculate cumulative viewshed from several viewsheds"""
gs.run_command("g.region", raster=maps, save="zoom_region", grow=250)
## faster shortcut, without reclass and HA_FILTER_NOISE:
# gs.run_command("r.series", input=maps, output="cumulative_viewshed", method="maximum")
gs.run_command("r.series", input=maps, output="cumulative_viewshed", method="sum")
gs.run_command("r.null", map="cumulative_viewshed", setnull=0)
# gs.run_command("r.colors", map="cumulative_viewshed", color="plasma")
gs.run_command("g.region", raster=maps, save="zoom_region", grow=250)
gs.mapcalc('cumulative_viewshed_bin = if(isnull(cumulative_viewshed), cumulative_viewshed, 1)')
gs.run_command(
"r.reclass.area", input="cumulative_viewshed_bin",
output="cumulative_viewshed_cleaned", mode="greater",
value=HA_FILTER_NOISE)
return "cumulative_viewshed_cleaned"
def save_viewshed(input_name, output_name):
"""Save viewshed to GeoTiff"""
gs.run_command("g.region", region="zoom_region")
gs.run_command(
"r.out.gdal",
flags="fmt",
input=input_name,
# input="cumulative_viewshed",
output=f"{PROJECT_FOLDER}out/{output_name}.tif",
format="GTiff",
overviews="4")
def cleanup_grass_project():
"""Clean up grass project layers (regions, raster, vector)"""
gs.run_command("g.remove",
type="vector",
pattern="pot_*",
flags="f")
gs.run_command("g.remove",
type="region",
pattern="zoom_*",
flags="f")
gs.run_command("g.remove",
type="raster",
pattern="viewshed_*",
flags="f")
gs.run_command("g.remove",
type="raster",
pattern="cumulative_*",
flags="f")
def remove_stopareas(input_name: str):
"""Remove stopareas from map"""
gs.mapcalc(
f'cumulative_viewshed_stopareas = if(isnull(stop_areas_raster), {input_name}, 0)')
gs.run_command("r.null", map="cumulative_viewshed_stopareas", setnull=0)
return "cumulative_viewshed_stopareas"
def display_raster_map(raster_layer):
"""Display raster map in jupyter as static image"""
viewshed_map = gj.Map(saved_region="zoom_region")
viewshed_map.d_rast(map="dsm")
viewshed_map.d_rast(map=raster_layer)
display(viewshed_map.show())
cleanup_grass_project()
%%bash
g.list type=raster,vector,region -m -t
get_shapes_list = list(Path("POT_WEA").glob('*.shp'))
len(get_shapes_list)
Process all viewsheds for all shapefiles
%%time
cleanup_grass_project()
# create progress bar for list of shapefiles
it = tqdm(get_shapes_list, unit='shapefile')
for shapefile in it:
loop_time = time.time()
cleanup_grass_project()
# gs.run_command("g.region", raster="dsm", save="zoom_region", grow=250)
gs.run_command("g.region", raster="dsm")
print(f"Processing {shapefile.name}")
viewpoints = get_viewpoints(shapefile.name)
maps = create_viewshed_group(viewpoints)
gs.run_command("g.region", raster=maps, save="zoom_region")
clear_output(wait=True)
# display_raster_map(raster_layer="viewshed_1")
output_name = shapefile.name[:-4]
input_name = cumulative_viewshed(maps=maps, output_name=output_name)
layer_stopareas = remove_stopareas(input_name=input_name)
# display_raster_map(raster_layer=layer_stopareas)
save_viewshed(layer_stopareas, output_name)
cleanup_grass_project()
it.display()
zip results¶
Make sure that 7z is available (apt-get install p7zip-full)
!cd /home/jovyan/work/viewshed_analysis/grassdata/viewshed && \
7z a -tzip -mx=9 out/$(date +"%Y-%m-%d")_viewsheds.zip \
out/* -x!*.zip -x!out/.ipynb_checkpoints \
-y > /dev/null
Create Notebook HTML¶
os.chdir("/home/jovyan/work/viewshed_analysis/notebooks")
!jupyter nbconvert --to html_toc \
--output-dir=../resources/html/ ./01_viewshedanalysis.ipynb \
--template=../nbconvert.tpl \
--ExtractOutputPreprocessor.enabled=False >&- 2>&- # create single output file