GIS/MEA582: Geospatial Modeling and Analysis

This is an unmaintained course material, please see current material at:

Viewshed and solar energy potential analysis

Resources:

Text files with site locations:

viewshed_points.txt

Start GRASS GIS

In 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) and for Accessible mapset create a new mapset (called e.g. HW_viewshed_solar) and click Start GRASS.

grass70

Change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI command 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 new 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.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.rast tower_los_lu30m
d.vect streets_wake
d.legend tower_los_lu30m at=25,50,1,3
d.out.file towerviewlu

We can also do visibility from former RedHat headquarters:

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.rast redhat_25_los
d.vect streets_wake
d.vect viewpoints size=10 col=red icon=basic/marker
d.out.file RHview

Use mapalgebra to compute landuse in the view, 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 RBC tower and RH headquarters.

Solar radiation 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.rast elevfacility_1m

Prepare input maps (slope and aspect):

r.slope.aspect 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 elevfacility_1m year=2001 month=12 day=22 hour=14 minute=25 sec=0 timezone=-5

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?

r.sun 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 compute and extract shadow area for 7.5 hr:

r.mapcalc "shadow_1m = if(isnull(incid_elevfac_1m), 1, null())"
r.colors shadow_1m color=grey
d.rast elevfacility_1m
d.rast shadow_1m
r.sun 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())"
d.rast shadow7_1m
d.out.file shadows

Solar radiation

Compute 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 view to see relation between terrain geometry and solar radiation.

r.sun 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 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.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,8867
d.out.file solar_summer

Calculate direct solar radiation 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 elev_ned_30m aspect=asp_ned_30m slope=slp_ned_30m
r.sun 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

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

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