Practical 7 Code Solution


Understanding GIS: Practical 7
@author jonnyhuck

Calculate a Viewshed using Bresenham's Line and Midpoint Circle Algorithm


New Topics:

from sys import exit
from geopandas import GeoSeries
from shapely.geometry import Point
from math import hypot, floor, ceil
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
from numpy import zeros, column_stack
from rasterio import open as rio_open
from matplotlib.colors import Normalize
from import ScalarMappable
from rasterio.plot import show as rio_show
from skimage.draw import line, circle_perimeter
from matplotlib.pyplot import subplots, savefig
from matplotlib_scalebar.scalebar import ScaleBar
from matplotlib.colors import LinearSegmentedColormap

def adjust_height(height, distance, earth_diameter=12740000, refraction_coefficient=0.13):
	* Adjust the apparant height of an object at a certain distance, accounting for the
	* 	curvature of the earth and atmospheric refraction
	* This is the Arc GIS version of the calculation for atmospheric refraction:
	* 	a = distance**2 / earth_diameter
	* 	return height - a + refraction_coefficient * a
	# this is the QGIS version of the calculation for atmospheric refraction
	return height - (distance**2 / earth_diameter) * (1 - refraction_coefficient)

def line_of_sight(r0, c0, height0, r1, c1, height1, radius, dem_data, d, output):
	 * Use Bresenham's Line algorithm to calculate a line of sight from one point to another point, 
	 *	returning a list of visible cells
	# init variable for biggest dydx so far (starts at -infinity)
	max_dydx = -float('inf')

	# loop along the pixels in the line (exclusing the first one)
	for r, c in column_stack(line(r0, c0, r1, c1))[1:]:

		# calculate distance travelled so far (in pixels)
		dx = hypot(c0 - c, r0 - r)
		# if we go too far, or go off the edge of the data, stop looping
		if dx > radius or not 0 <= r < d.height or not 0 <= c < d.width:

		# calculate the current value for dy / dx
		base_dydx = (adjust_height(dem_data[(r, c)], dx * d.res[0]) - height0) / dx
		tip_dydx = (adjust_height(dem_data[(r, c)] + height1, dx * d.res[0]) - height0) / dx
		# if the tip dydx is bigger than the previous max, it is visible
		if tip_dydx > max_dydx:
			output[(r, c)] = 1

		# if the base dydx is bigger than the previous max, update
		max_dydx = max(max_dydx, base_dydx)

	# return updated output surface
	return output

def viewshed(x0, y0, radius_m, observer_height, target_height, dem_data, d):
	* Use Midpoint Circle algorithm to determine endpoints for viewshed
	# convert origin coordinates to image space
	r0, c0 = d.index(x0, y0)

	# make sure that we are within the dataset
	if not 0 <= r0 < d.height or not 0 <= c0 < d.width:
		print(f"Sorry, {x0, y0} is not within the elevation dataset.")

	# convert the radius (m) to pixels
	radius_px = int(radius_m / d.res[0])

	# get the observer height (above sea level)
	height0 = dem_data[(r0, c0)] + observer_height
	# create output array at the same dimensions as data for viewshed
	output = zeros(dem_data.shape)

	# set the start location as visible automatically
	output[(r0, c0)] = 1

	# get pixels in the perimeter of the viewshed
	for r, c in column_stack(circle_perimeter(r0, c0, radius_px*3)):

		# calculate line of sight to each pixel, pass output and get a new one back each time
		output = line_of_sight(r0, c0, height0, r, c, target_height, radius_px, dem_data, d, output)

	# return the resulting viewshed
	return output

# open the elevation data file
with rio_open("../data/helvellyn/Helvellyn-50.tif") as d:

	# read the data out of band 1 in the dataset
	dem_data =

	# set origin for viewshed
	x0, y0 = 332000, 514000

	# calculate the viewshed
	output = viewshed(x0, y0, 5000, 1.8, 100, dem_data, d)

	# output image
	fig, my_ax = subplots(1, 1, figsize=(16, 10))
	my_ax.set_title("Viewshed Analysis")

	# draw dem
		transform = d.transform,
		cmap = 'viridis',

	# draw dem as contours
		transform = d.transform,
		colors = ['white'],
		linewidths = [0.5],

	# add viewshed
		cmap = LinearSegmentedColormap.from_list('binary_viewshed', [(0, 0, 0, 0), (1, 0, 0, 0.5)], N=2)

	# add origin point
	GeoSeries(Point(x0, y0)).plot(
		ax = my_ax,
		markersize = 60,
		color = 'black',
		edgecolor = 'white'

	# add a colour bar
	fig.colorbar(ScalarMappable(norm=Normalize(vmin=floor(dem_data.min()), vmax=ceil(dem_data.max())), cmap='viridis'), ax=my_ax, pad=0.01)

	# add north arrow
	x, y, arrow_length = 0.97, 0.99, 0.1
	my_ax.annotate('N', xy=(x, y), xytext=(x, y-arrow_length),
		arrowprops=dict(facecolor='black', width=5, headwidth=15),
		ha='center', va='center', fontsize=20, xycoords=my_ax.transAxes)

	# add scalebar
	my_ax.add_artist(ScaleBar(dx=1, units="m", location="lower right"))

	# add legend for point
			Patch(facecolor=(1, 0, 0, 0.5), edgecolor=None, label=f'Visible Area'),
			Line2D([0], [0], marker='o', color=(1,1,1,0), label='Viewshed Origin', markerfacecolor='black', markersize=8)
		], loc='lower left')

	# save the result
	savefig('./out/7.png', bbox_inches='tight')