<jupyter_start><jupyter_text>Viewshed and solar energy potential analysisResources:* [GRASS GIS overview and manual](http://grass.osgeo.org/grass72/manuals/index.html)* [ GRASSbook](http://www.grassbook.org/)Text files with site locations:* [viewshed_points.txt](data/viewshed_points.txt) Start GRASS GISIn startup pannel set GIS Data Directory to path to datasets,for example on MS Windows, `C:\Users\myname\grassdata`.For Project location select nc_spm_08_grass7 (North Carolina, State Plane, meters) andfor Accessible mapset create a new mapset (called e.g. HW_viewshed_solar) andclick Start GRASS.<jupyter_code># This is a quick introduction into Jupyter Notebook.
# Python code can be excecuted like this:
a = 6
b = 7
c = a * b
print "Answer is", c
# Python code can be mixed with command line code (Bash).
# It is enough just to prefix the command line with an exclamation mark:
!echo "Answer is $c"
# Use Shift+Enter to execute this cell. The result is below.
import os
import sys
import subprocess
from IPython.display import Image

# create GRASS GIS runtime environment
gisbase = subprocess.check_output(["grass", "--config", "path"]).strip()
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))

# do GRASS GIS imports
import grass.script as gs
import grass.script.setup as gsetup

# set GRASS GIS session data
rcfile = gsetup.init(gisbase, "/home/jovyan/grassdata", "nc_spm_08_grass7", "user1")
# default font displays
os.environ['GRASS_FONT'] = 'sans'
# overwrite existing maps
os.environ['GRASS_OVERWRITE'] = '1'
gs.set_raise_on_error(True)
gs.set_capture_stderr(True)
# set display modules to render into a file (named map.png by default)
os.environ['GRASS_RENDER_IMMEDIATE'] = 'cairo'
os.environ['GRASS_RENDER_FILE_READ'] = 'TRUE'
os.environ['GRASS_LEGEND_FILE'] = 'legend.txt'<jupyter_output><empty_output><jupyter_text>Change working directory:_Settings_ > _GRASS working environment_ > _Change working directory_ > select/create any directoryor type `cd` (stands for change directory) into the GUI_Console_ and hit Enter:<jupyter_code># a proper directory is already set, download files
import urllib
urllib.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/viewshed_points.txt", "viewshed_points.txt")<jupyter_output><empty_output><jupyter_text>Download all text files (see above)to the selected directory. Now you can use the commands from the assignment requiring the text filewithout the need to specify the full path to the file. Viewshed analysisCompute viewshed from a new 32 story tower located in downtown.<jupyter_code>gs.parse_command('g.region', raster="elev_ned_30m", flags='apg')
gs.run_command('v.in.ascii', input="viewshed_points.txt", output="viewpoints", separator=",", z="3", flags='z')
gs.run_command('r.viewshed', input="elev_ned_30m", output="tower_165_los", coordinates="642212,224767", observer_elevation="165", max_distance="10000")<jupyter_output><empty_output><jupyter_text>Display result on shaded relief:<jupyter_code>gs.run_command('d.his', hue="tower_165_los", intensity="elevation_shade")
gs.run_command('d.vect', map="streets_wake")
gs.run_command('d.vect', map="viewpoints", size="10", color="orange", icon="basic/marker")
gs.run_command('d.barscale', at="3,6")
gs.run_command('d.legend', raster="tower_165_los", at="50,90,2,6")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Find out what is the landuse within the view using map algebra:<jupyter_code>gs.mapcalc("tower_los_lu30m = if(tower_165_los,landclass96,null())")
gs.run_command('r.colors', map="tower_los_lu30m", raster="landclass96")
print gs.read_command('r.category', map="tower_los_lu30m", raster="landclass96")
print gs.read_command('r.report', map="tower_los_lu30m", unit="p,h", flags='n')<jupyter_output><empty_output><jupyter_text>Display only the following layers and save result:<jupyter_code>gs.run_command('d.rast', map="tower_los_lu30m")
gs.run_command('d.vect', map="streets_wake")
gs.run_command('d.legend', raster="tower_los_lu30m", at="25,50,1,3")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>We can also do visibility from former RedHat headquarters:<jupyter_code>gs.run_command('r.viewshed', input="elev_ned_30m", output="redhat_25_los", coordinates="638898,224528", observer_elevation="25", max_distance="10000")<jupyter_output><empty_output><jupyter_text>Display only the following layers and save result:<jupyter_code>gs.run_command('d.rast', map="redhat_25_los")
gs.run_command('d.vect', map="streets_wake")
gs.run_command('d.vect', map="viewpoints", size="10", col="red", icon="basic/marker")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Use mapalgebra to compute landuse in the view,assign the visible land use map land use colors and category labelsusing _r.colors_ and _r.category_ (see previous example)and use _r.report_ -n to compare the size and land use composition withinviewshed from the RBC tower and RH headquarters. Solar radiation analysisSet the region and add the planned building to the DEM, we will use this new DEM for the analyses.Remove all layers and zoom to the region.<jupyter_code>gs.parse_command('g.region', region="rural_1m", flags='pg')
gs.mapcalc("elevfacility_1m = if(isnull(facility), elev_lid792_1m,138.)")
gs.run_command('r.colors', map="elevfacility_1m", color="elevation")
gs.run_command('d.erase')
gs.run_command('d.rast', map="elevfacility_1m")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Prepare input maps (slope and aspect):<jupyter_code>gs.run_command('r.slope.aspect', elevation="elevfacility_1m", aspect="aspect_elevfac_1m", slope="slope_elevfac_1m")<jupyter_output><empty_output><jupyter_text>Incidence angles and cast shadowsCompute the sun position on Dec. 22 at 2:25pm, EST (no map output expected):<jupyter_code>gs.run_command('r.sunmask', elevation="elevfacility_1m", year="2001", month="12", day="22", hour="14", minute="25", sec="0", timezone="-5", flags='s')<jupyter_output><empty_output><jupyter_text>Calculate incidence angles including cast shadows.We assign histogram equalized color table - can you explain why?(hint: try the same color table without -e).What is the value on the roof? How is it related to day/time?<jupyter_code>gs.run_command('r.sun', elevation="elevfacility_1m", aspect="aspect_elevfac_1m", slope="slope_elevfac_1m", incidout="incid_elevfac_1m", day="356", time="14.416667")
gs.parse_command('r.info', map="incid_elevfac_1m", flags='g')
gs.run_command('r.colors', map="incid_elevfac_1m", co="bcyr", flags='e')
gs.run_command('d.rast', map="incid_elevfac_1m")
gs.run_command('d.legend', raster="incid_elevfac_1m", at="25,50,1,3")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Extract the cast shadow area for 14.4 hr andcompute and extract shadow area for 7.5 hr:<jupyter_code>gs.mapcalc("shadow_1m = if(isnull(incid_elevfac_1m), 1, null())")
gs.run_command('r.colors', map="shadow_1m", color="grey")
gs.run_command('d.rast', map="elevfacility_1m")
gs.run_command('d.rast', map="shadow_1m")
gs.run_command('r.sun', elevation="elevfacility_1m", aspect="aspect_elevfac_1m", slope="slope_elevfac_1m", incidout="incid_elevfac7_1m", day="356", time="7.50")
gs.mapcalc("shadow7_1m = if(isnull(incid_elevfac7_1m), 1, null())")
gs.run_command('d.rast', map="shadow7_1m")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Solar radiationCompute global (beam+diffuse+refl) radiation for entire day during summer and winter solstice.Display the radiation maps and also insolation time maps insol_time using same commands.Optionally display the radiation maps draped over elevation elevfacility_1m in 3D viewto see relation between terrain geometry and solar radiation.<jupyter_code>gs.run_command('r.sun', elevation="elevfacility_1m", aspect="aspect_elevfac_1m", slope="slope_elevfac_1m", day="356", glob_rad="g356", insol_time="its356")
gs.run_command('r.colors', map="g356", co="gyr", flags='e')
gs.run_command('r.sun', elevation="elevfacility_1m", aspect="aspect_elevfac_1m", slope="slope_elevfac_1m", day="172", glob_rad="g172", insol_time="its172")
gs.run_command('r.colors', map="g172", co="gyr", flags='e')
gs.run_command('d.rast', map="g356")
gs.run_command('d.legend', raster="g356", at="25,50,1,3")
Image(filename="map.png")
gs.run_command('d.rast', map="g172")
gs.run_command('d.legend', raster="g172", at="25,50,1,3", range="8800,8867")
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Calculate direct solar radiation and insolation time for a larger region.Try to find good color tables (custom, hist. equalized, to see the pattern).<jupyter_code>gs.parse_command('g.region', raster="elev_ned_30m", flags='pg')
gs.run_command('r.slope.aspect', elevation="elev_ned_30m", aspect="asp_ned_30m", slope="slp_ned_30m")
gs.run_command('r.sun', elevation="elev_ned_30m", aspect="asp_ned_30m", slope="slp_ned_30m", linke_value="2.5", albedo_value="0.2", beam_rad="b356", diff_rad="d356", refl_rad="r356", insol_time="it356", day="356")<jupyter_output><empty_output><jupyter_text>Zoom to the new computational region and display the following layers:<jupyter_code>gs.run_command('d.rast', map="b356")
gs.run_command('d.legend', raster="b356", at="2,30,2,6")
Image(filename="map.png")
gs.run_command('d.rast', map="it356")
gs.run_command('d.legend', raster="it356", at="2,30,2,6")
Image(filename="map.png")
# end the GRASS session
os.remove(rcfile)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14655172413793102, 'fraction_numerical': 0.03353658536585366, 'mean_word_length': 8.726993865030675, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 68, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '74187', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3586, 'n_tokens_neox': 3370, 'n_words': 725, 'path': 'notebooks/viewshed_solar_python.ipynb', 'repo_name': 'Geethen/geospatial-modeling-course-jupyter', 'size': '9512'}