library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(readxl)
library(RColorBrewer)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(viridis)
## Loading required package: viridisLite
library(tmap)
library(here)
## Warning: package 'here' was built under R version 4.2.2
## here() starts at C:/temp/class2122/CP8853
source("functions.R")
# climate data is in wgs84 (EPSG 4326)
#    other data will need to be re-projected
vwgs84 = 4326
rwgs84 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
usst1 = st_read("cb_2018_us_state_20m.shp") %>%
  tolow() %>%
  # restrict analysis to coterminous US
  filter(name != 'Puerto Rico' &
           name != 'American Samoa' &
           name != 'United States Virgin Islands' &
           name != 'Commonwealth of the Northern Mariana Islands' &
           name != 'Alaska' &
           name != 'Hawaii' &
           name != 'Guam') %>% 
  st_transform(vwgs84)
## Reading layer `cb_2018_us_state_20m' from data source 
##   `C:\temp\class2122\CP8853\cb_2018_us_state_20m.shp' using driver `ESRI Shapefile'
## Simple feature collection with 52 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
## Geodetic CRS:  NAD83
usco1 = st_read("cb_2018_us_county_20m.shp") %>%
  tolow()  %>% 
  mutate(geoid = paste0("g",geoid)) %>% 
  st_transform(wgs84) %>% 
  # restrict analysis to coterminous US
  st_intersection(usst1)
## Reading layer `cb_2018_us_county_20m' from data source 
##   `C:\temp\class2122\CP8853\cb_2018_us_county_20m.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3220 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
## Geodetic CRS:  NAD83
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
roads = st_read("primaryroads4.shp") %>%
  tolow() %>% 
  # limit to Interstates
  filter(rttyp == "I") %>% 
  st_transform(vwgs84)
## Reading layer `primaryroads4' from data source 
##   `C:\temp\class2122\CP8853\primaryroads4.shp' using driver `ESRI Shapefile'
## Simple feature collection with 12105 features and 6 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -123.3938 ymin: 25.74903 xmax: -67.78126 ymax: 49.00236
## Geodetic CRS:  NAD83
gast1 = usst1 %>% 
  filter(statefp=='13') 

gaco1 = usco1 %>% 
  filter(statefp == "13") 

fuldekco1 = gaco1 %>% filter(name %in% c("Fulton", "Dekalb"))
taylorco1 = gaco1 %>% filter(name %in% c("Taylor"))

#Task#2- Drawing tmap

tm_shape(gaco1) +
  tm_polygons(col="pink") +
tm_shape(gast1) +
  tm_borders(lwd=3)

#Task#3- Map of GA using a color ramp

tm_shape(gaco1) +
  tm_polygons(col="aland",palette="PuRd") +
tm_shape(gast1) +
  tm_borders(lwd=3)

#Task#4- Reading geotif raster files #Explanation from Professor Drummond’s Lab: #GHI stands for global horizontal irradiance. It is the basic measure of solar radiation used for flat-plane solar collectors such as solar panels. The GHI calculation includes factors, such as clound cover, that reduce the amount of effective solar radition.

#GHI maps usually show kWh per square meter per day (kWh/m2/d). The GHI is usually a mean monthly or annual value, and is often averaged over multiple years.

#The raster package command raster() will read various types of raster files. It will use the three letter file extension to determine the type of the input file.

#The following commands read national U.S. GHI geotif file developed by the National Renewables Energy Laboratory (NREL) and set the CRS to WGS84. Note: the raster package uses a proj4 CRS formulation to set CRS information with either

#an “+init=” value for an EPSG designation or #a full set of proj4 “+” parameters, and #unlike sf it can not accept a WKT CRS format

usghi = raster("nsrdb3_ghi.tif")
crs(usghi) = rwgs84

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(usghi) +
  tm_raster()
## stars object downsampled to 1395 by 717 cells. See tm_shape manual (argument raster.downsample)

#Obtaining raster basic info

usghi
## class      : RasterLayer 
## dimensions : 2025, 3940, 7978500  (nrow, ncol, ncell)
## resolution : 0.04, 0.04  (x, y)
## extent     : -180, -22.4, -21.01, 59.99  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : nsrdb3_ghi.tif 
## names      : nsrdb3_ghi

#Task#5- Creating a sub-raster #To generate a local raster file that includes only Georgia we mustapply a vector polygon mask to select only raster cells inside the Georgia polygon and #reset the raster size by cropping the bounding box to the Georgia polygon bounding box.

#To save a bit of coding work we can apply the raster package’s mask and crop functions in the same step by nesting the mask function inside the crop function.

source("functions.R")
gaghi = crop(mask(usghi,gast1),gast1)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(gaghi) +
  tm_raster() +
  tm_layout(legend.position=c("right","top"))

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(gaghi) +
  tm_raster(alpha=.8, pal="inferno",
            title="GHI: kWh/m2/d") +
tm_shape(gaco1) +
  tm_borders() +  
tm_layout(legend.position=c("right", "top"),
          legend.bg.color="white")

#Task#6- Building an interactive map

tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap") +
tm_shape(gaghi) +
  tm_raster(alpha=.8, pal="inferno",
            title="GHI: kWh/m2/d") +
tm_shape(gaco1) +
  tm_borders() +
tm_shape(roads) +
  tm_lines(col="cyan")

#Interactive Map Using Esri.DeLorme basemap

tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("Esri.DeLorme") +
tm_shape(gaghi) +
  tm_raster(alpha=.8, pal="inferno",
            title="GHI: kWh/m2/d") +
tm_shape(gaco1) +
  tm_borders()

#Task#7- Averaging raster values for polygons

mapvalue1 = raster::extract(gaghi,gaco1,fun=mean,na.rm=TRUE)
mapvalue2 = as.data.frame(mapvalue1)
gaco2 = bind_cols(gaco1,mapvalue2) %>% 
  rename(ghi=V1) %>% 
  dplyr::select(name, ghi)
mapvalue1 = raster::extract(gaghi,gaco1,fun=mean,na.rm=TRUE)
mapvalue2 = as.data.frame(mapvalue1)
gaco2 = bind_cols(gaco1,mapvalue2) %>% 
  rename(ghi=V1) %>% 
  dplyr::select(name, ghi)
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(gaco2) +
  tm_polygons("ghi",palette="PuRd")

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(gaco2) +
  tm_polygons("ghi",palette="PuRd")

#Task#8- Raster statistics

summary(gaghi)
##         nsrdb3_ghi
## Min.         4.416
## 1st Qu.      4.752
## Median       4.848
## 3rd Qu.      4.920
## Max.         5.088
## NA's      4628.000
tmap_mode("plot")
## tmap mode set to plotting
hist(gaghi)

boxplot(gaghi)

#Task#9- Advanced Terrain functions

tmap_mode("plot")
## tmap mode set to plotting
e1 = raster("gadem2.tif")

tm_shape(e1) +
 tm_raster() +
tm_shape(gaco1) +
 tm_borders() +
tm_layout(legend.position=c("right","top"),
          legend.bg.color="white")
## stars object downsampled to 948 by 1056 cells. See tm_shape manual (argument raster.downsample)

#Question 1: How do you get the CRS information for the raster?

#Question 2: What are the xy units?

#Question 3: What are the z units? (You may need to do some quick Internet research to answer this question.)

#Question 4: If the vertical units are meters and the horizontal ones are feet, what could be the problem in calculating slope?

#Task#10- Creating slope and aspect rasters #We can now use the terrain function to create slope and aspect maps. Aspect is the direction that sloped land faces.

#The tmap commands use the fuldekco1 file to zoom to the Atlanta area. Note Stone Mountain to the east, and the Chattahoochee River flowing to the southwest along Fulton County’s eastern border.

s1 = terrain(e1, opt="slope", unit="degrees", neighbors=4)
a1 = terrain(e1, opt="aspect", unit="degrees", neighbors=4)

tm_shape(fuldekco1) +
  tm_borders() +
tm_shape(s1) +
 tm_raster(palette="YlGn",breaks=c(0,1,3,5,10,15,40)) +
tm_shape(gaco1) +
 tm_borders() +
tm_layout(legend.position=c("left","top"),
          legend.bg.color="white")
## stars object downsampled to 948 by 1056 cells. See tm_shape manual (argument raster.downsample)
## Warning: Values have found that are higher than the highest break