Assignment 02

Does the sun shine there?

Deadline is 2022-12-08 17:00

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

This is an individual assignment.

• The DSM to use: AHN3
• If you choose Python
• If you choose C++
• Good to know
• Marking
• What to submit and how to submit it

This assignment is inspired from the ShadowAnimation module of WhiteboxTools, which you can see above.

You have to implement one part of this algorithm: the line-of-sight between a location around Delft and the sun. Basically we want to know whether at a given location at a given time the sun shines (on the ground), or not. Notice that we assume Delft to be cloud- and rain-free 😎

You have to implement a variation of the line-of-sight algorithm, as explained in the book (Chapter 9).

The DSM to use: AHN3

Your code should work for any DSM 0.5m in the AHN3.

Since each tile is pretty big (roughly 500MB) I cropped from tile 37EN2 a subset around BK-City, use it to test your implementation, it’s in the GitLab repository in 2 formats (hw/02/ahn3_data/ahn3_dsm50cm_bk_small.tif/asc).

You’ll notice that the DSM has several no_data values, there is no need to fill them, just assume that there is nothing there.

If you choose Python

Code is in the /hw/02/python/ folder of the GitLab repository of the course.

We give you the skeleton of the Python program:

1. geo1015_hw02.py is the main(). You are not allowed to modify this file! This main calls the is_sunny() function in mycode_hw02.py.
2. mycode_hw02.py is where all your code should go. The is_sunny() function must be completed, and you are not allowed to change the input parameters. You are of course allowed to add any other functions you want in that unit (please do not include any other unit you code, add your code to that single file). You are only allowed to import modules from the Python standard library (anything you installed through pip is thus not allowed, except the ones already defined there).

You have to write the code for the is_sunny() function in mycode_hw02.py.

It returns True if p (the location) at date-time dt (Delft local time, not UTC) has the ground illuminated by the sun, False otherwise. The ground is illuminated by the sun if the line-of-sight between that point and the sun is not obstructed by any obstacles (buildings or lampposts or trees or anything in the AHN3 DSM dataset).

However, if p is outside the spatial extent of the input raster (given in EPSG:28992), then the function needs to throw an Exception: raise Exception("Point given is outside the extent of the dataset").

Also, if p is located where there is no elevation value (no_data), then the function needs to throw an Exception: raise Exception("Point given has no_data value").

You’ll need to install the following libraries, which are not part of the standard library of Python:

1. suncalc
2. pytz
3. pyproj

What is necessary to install is thus:

pip install suncalc 
pip install pytz
pip install pyproj
pip install rasterio
pip install numpy

To run the code:

python geo1015_hw02.py ../../ahn3_data/ahn3_dsm50cm_bk_small.tif 85000 450000 '2022-11-30 11:46'

If you choose C++

Code is in the /hw/02/cpp/ folder of the GitLab repository of the course.

We give you the skeleton of the C++ program:

1. /src/main.cpp has the main() function. The function is_sunny() is in this unit, so you need to complete it. You are of course allowed to add any other functions you want.
2. /src/DatasetASC.h & /src/DatasetASC.cpp is the class where the input .asc file is read and stored.
3. ../ahn3_data/: contains 1 example of a the AHN3 DSM 0.5m raster in .asc format. (The .asc format is used because reading GeoTIFF with C++ could require installing complex libraries and under Windows it could be difficult)

You have to write the code for the is_sunny() function in main.cpp.

Since suncalc and the reprojection of CRS (with proj) require rather complex libraries to install, I took the same Python code that I give to students who chose Python, and the Python unit is called from C++ and the return value read (see the example in the code) (it’s not pretty but it works…).

DatasetASC.h has a few functions to allow you to perform the same as for those using rasterio, eg the file xy2rc() and rc2xy() to convert from row/col to x/y and vice-versa.

You don’t need to install anything to compile and run the code (with the console):

mkdir build
cd build
cmake ..
make
./hw02 ../../ahn3_data/ahn3_dsm50cm_bk_small.asc 86000 445000 '2022-09-26 11:35'

At run-time, the code calls the Python code /python_bin/sunpos.py to run this you need to install those:

pip install suncalc 
pip install pytz
pip install pyproj

If you want to use CLion, just open the cpp folder as a project and all will be fine. To run the code with the 4 arguments, you do this (notice the double-quotes for the date-time):

Good to know

• suncalc uses only UTC time (Coordinated Universal Time), both as input and output. Delft is not UTC so you need to convert back and forth
• suncalc uses only EPSG:4326 and the coordinates used as input have to be in EPSG:28992, so you need to convert them
• the time it takes for your program to terminate has no influence on the marks
• the marking will be automatic, I will test your mycode_hw02.py/main.cpp with different inputs in Delft (so the time zone is always that of Delft) and verify whether you return the correct values
• 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 here, and there. The second method is actually the same as the conversion vector-raster as seen in Lecture 10 of GEO1002, and it is already implemented in rasterio. We give you an example in the Python code (see its docs)

Marking

What to submit and how to submit it

You have to submit one file: my_code_hw02.py (or main.cpp) where your name is clearly identified at the top.

Oh, and no report is necessary.

Upload the file file to this Dropbox page.

[last updated: 2022-11-30 12:05]