GIS/MEA582: Geospatial Modeling and Analysis

Viewshed and solar energy potential analysis

Resources:

Text files with site locations:

viewshed_points.txt

Start GRASS GIS

Create a new mapset (called e.g. HW_viewshed_solar) in nc_spm_08_grass7 project and change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter:

cd

Download all text files (see above) to the selected directory. Now you can use the commands from the assignment requiring the text file without the need to specify the full path to the file.

Viewshed analysis

Compute viewshed from a 32 story tower located in downtown.

g.region raster=elev_ned_30m -ap
v.in.ascii -z input=viewshed_points.txt output=viewpoints separator=, z=3
r.viewshed elev_ned_30m output=tower_165_los coordinates=642212,224767 observer_elevation=165 max_distance=10000

Display result on shaded relief:

d.erase
d.his hue=tower_165_los intensity=elevation_shade
d.vect streets_wake
d.vect viewpoints size=10 color=orange icon=basic/marker
d.barscale at=3,6
d.legend tower_165_los at=50,90,2,6
d.out.file towerview_angle

Find out what is the landuse within the view using map algebra:

r.mapcalc "tower_los_lu30m = if(tower_165_los,landclass96,null())"
r.colors tower_los_lu30m raster=landclass96
r.category tower_los_lu30m raster=landclass96
r.report -n tower_los_lu30m unit=p,h

Display only the following layers and save result:

d.erase
d.rast tower_los_lu30m
d.vect streets_wake
d.legend -c -b tower_los_lu30m at=25,50,1,3
d.out.file towerviewlu

We can also map visibility from former RedHat headquarters on Centennial Campus:

r.viewshed elev_ned_30m output=redhat_25_los coordinates=638898,224528 observer_elevation=25 max_distance=10000

Display only the following layers and save result:

d.erase
d.rast redhat_25_los
d.vect streets_wake
d.vect viewpoints size=10 col=red icon=basic/marker
d.out.file RHview

Use map algebra to compute landuse in within the viewshed, assign the visible land use map land use colors and category labels using r.colors and r.category (see previous example) and use r.report -n to compare the size and land use composition within viewshed from the downtown tower and the building on Cenetennial Campus that served as RedHat headquarters.

Solar irradiation analysis

Set 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.

g.region rural_1m -p
r.mapcalc "elevfacility_1m = if(isnull(facility), elev_lid792_1m,138.)"
r.colors elevfacility_1m color=elevation
d.erase
d.rast elevfacility_1m

Prepare input maps (slope and aspect):

r.slope.aspect elevation=elevfacility_1m aspect=aspect_elevfac_1m slope=slope_elevfac_1m

Incidence angles and cast shadows

Compute the sun position on Dec. 22 at 2:25pm, EST (no map output expected):

r.sunmask -s elevation=elevfacility_1m year=2001 month=12 day=22 hour=14 minute=25 sec=0 timezone=-5

Calculate incidence angles including cast shadows on Dec. 22 at 2:25pm, EST.
We assign histogram equalized color table - can you explain why?
What is the value of the incidence angle on the roof? How is it related to solar altitude? (hint: see the output of r.info and query the incidence angle map on the building)

r.sun elevation=elevfacility_1m aspect=aspect_elevfac_1m slope=slope_elevfac_1m incidout=incid_elevfac_1m day=356 time=14.416667
r.info incid_elevfac_1m
r.colors -e incid_elevfac_1m co=bcyr
d.rast incid_elevfac_1m
d.legend incid_elevfac_1m at=25,50,1,3
d.out.file incid_angle

Extract the cast shadow area for 14.4 hr and then compute an incidence angle map and extract shadow area for the morning at 7.5 hr same day:

r.mapcalc "shadow_1m = if(isnull(incid_elevfac_1m), 1, null())"
r.colors shadow_1m color=grey
d.erase
d.rast elevfacility_1m
d.rast shadow_1m
d.out.file shadows14hr
r.sun elevation=elevfacility_1m aspect=aspect_elevfac_1m slope=slope_elevfac_1m incidout=incid_elevfac7_1m day=356 time=7.50
r.mapcalc "shadow7_1m = if(isnull(incid_elevfac7_1m), 1, null())"
r.colors shadow7_1m color=grey
d.erase
d.rast elevfacility_1m
d.rast shadow7_1m
d.out.file shadows7hr

How does the cast shadow area change between the morning and early afternoon on the day of winter solstice? (hint, use r.report) Is there any overlap?

Solar irradiation

Compute global (beam+diffuse+refl) irradiation for entire day during summer and winter solstice.
Display the irradiation maps and insolation time maps with legends to compare the range of values

r.sun elevation=elevfacility_1m aspect=aspect_elevfac_1m slope=slope_elevfac_1m day=356 glob_rad=g356 insol_time=its356
r.colors g356 co=gyr -e
r.sun elevation=elevfacility_1m aspect=aspect_elevfac_1m slope=slope_elevfac_1m day=172 glob_rad=g172 insol_time=its172
r.colors g172 co=gyr -e
d.erase
d.rast g356
d.legend g356 at=25,50,1,3
d.out.file solar_winter
d.rast g172
d.legend g172 at=25,50,1,3 range=8800,9000
d.out.file solar_summer
d.erase
d.rast its356
d.legend its356 at=25,50,1,3
d.out.file solartime_winter
d.rast its172
d.legend its172 at=25,50,1,3 
d.out.file solartime_summer

Discuss the difference between the daily solar irradiation values and insolation time during the days of summer and winter solstice. What are the units for solar irradiation? What are the implications for using solar irradiation for generating sufficient energy throughout the year?

Optionally display the irradiation maps draped over elevation elevfacility_1m in 3D view to see relation between terrain geometry and solar irradiation.

Optional: Calculate direct solar irradiation and insolation time for a larger region.

Try to find good color tables (custom, hist. equalized, to see the pattern).

g.region raster=elev_ned_30m -p
r.slope.aspect elevation=elev_ned_30m aspect=asp_ned_30m slope=slp_ned_30m
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
r.colors b356,d356,r356 co=bcyr -e

Zoom to the new computational region and display the following layers:

d.erase
d.rast b356
d.legend b356 at=2,30,2,6
d.out.file beamrad356
d.rast it356
d.legend it356 at=2,30,2,6
d.out.file insol_time