Homework 03: Computing a viewshed

 

Deadline is 2019-01-14 13:00.

Late submission? 10% will be removed for each day that you are late.

You’re allowed for this assignment to work in a group of 2 (and thus submit only one report for both of you). If you prefer to work alone it’s also fine (but we don’t recommend it).

---

• Overview
• What you are given to start
• The output format
• You need to design a strategy yourself, and document it in a report
• Good to know
• What to submit and how to submit it

---

Overview

The aim of this assignment is to design and implement viewshed computation for grids in Python. Eventually, you should be able to compute what parts of Rio de Janeiro can be seen by Cristo Redentor, the famous 38m high statue shown in the image above. Or actually, what he would see if he could turn around his head…

The handout provides more information on how to do viewshed computation. You will have to implement the algorithm so that several viewpoints are taken into account. The resulting viewshed is what can be seen from any of the viewpoints, ie if one DTM cell is visible only from one viewpoint, then it is visible.

What you are given to start

code_hw03.zip

We give you the skeleton of the Python program:

1. geo1015_hw03.py is the main(). You are not allowed to modify this file!
2. mycode_hw03.py is where all your code should go. The function output_viewshed() receives all the relevant information and writes to the disk the output viewshed, as defined below. You are of course allowed to add any other functions you want. You are only allowed to import modules from the Python standard library (anything you installed through pip is thus not allowed, except rasterio, scipy and numpy).
3. params.json: a very simple JSON file that defines what the input file is, the viewpoint(s), with what parameters the viewpoints will be computed, and where to store the output. The “height” is the height of the viewpoint above the surface (2m for a person thus).
4. /data/: contains 2 raster terrains. For /data/CristoRendetor/, the location of the viewpoint is that of Cristo Redentor himself.
5. example_out.tif: an example of a GeoTIFF output file that your program will need to create.

The input file params.json given doesn’t need to be validated, and the input/output and parameters are also valid.

The output format

The output file has to be a GeoTiff (.tif), with the same CRS/extent/resolution as the input, but with type “Byte - Eight bit unsigned integer” (the code given shows how to do this).

There are 4 possible values for a pixel:

• 0: not visible from the viewpoint(s) (but inside the max-distance/horizon zone(s))
• 1: visible from the viewpoint(s)
• 2: the pixel contains a viewpoint
• 3: the pixel is outside the max-distance/horizon zone(s)

You need to design a strategy yourself, and document it in a report

There are several ways to solve this problem, and you are free to design the solution you want. Doing it brute-force, as described in the handout, is fine.

You are allowed to make approximations to speed up the operations performed, but these should be documented in your report. For instance, a grid of ~30m resolution brings a certain uncertainty, so one is allowed to approximate some parts when extracting which pixels represent the line segments between the viewpoint and a target point.

Your report does not need to have an introduction. It should only describe your algorithm at a high level (using pseudo-code, or a flowchart). Start by presenting the big steps, and then describe briefly how you implemented them. If assumptions were made to speed up the code, describe and document in the report the consequences of these. The report should ideally be 4-5 pages.

Also, you should report how long it took to generate the 2 files (for params.json and params2.json). The time is printed automatically when the process finishes.

Good to know

• you can assume that the viewpoints are at the centre of the pixels in which they are located (thus so you can shift the viewpoints slightly)
• the CRS of the input grids will always be in meters
• the viewpoints can “look” in all directions
• the viewpoints are guaranteed to be inside the cells (they are given in xy real-coordinates)
• identifying all the grid cells intersecting a line will be slow if the intersection with each must be calculated. Some ideas on how to perform this are there, and there. The second method is actually the same as the conversion vector-raster as seen in Lecture 10 of GEO1002, and it is implemented in rasterio. We give you an example in the given code.
• the QGIS plugin Profile Tool is useful to debug and to assess whether you obtain sensible results

What to submit and how to submit it

You have to submit a total of 4 files:

1. the Python file my_code_hw03.py, where your names are clearly identified at the top of the file.
2. the output GeoTIFF file out-tasmania.tif (for params.json)
3. the output GeoTIFF file out-cristo.tif (for params2.json)
4. a report in PDF format (no Word file please)

Do not submit your assignment by email, but upload the requested files zipped to this Dropbox file request page. Make sure that you put the full name of one of the members of the team (only one is sufficient).

You’ll get a confirmation email when everything has been successfully uploaded, keep it in case things go wrong.

[last updated: 2018-12-19]