In this exercise we will plot the west wide sit suitability data using simple features and tmap packages.

# install sf package if you haven't already
# install.packages("sf")
# load library
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
# Check working directory to set this up correctly
# Commands  getwd(), my.dir <- getwd(), list.files(my.dir, recursive= TRUE) 
# getwd()
# [1] "C:/Users/Emily Leslie/Documents/earth-analytics/data"

# import data

# To import shapefiles we use the sf library function st_read(). st_read() requires the file path to the shapefile.
ct <- st_read("continent/continent.shp")
## Reading layer `continent' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\continent\continent.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -16900910 ymin: -6972041 xmax: 16900900 ymax: 15297990
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
# Reading layer `continent' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\continent\continent.shp' using driver `ESRI Shapefile'
# Simple feature collection with 8 features and 5 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -16900910 ymin: -6972041 xmax: 16900900 ymax: 15297990
# epsg (SRID):    NA
# proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

#install packages and load tmap library if you haven't already 
# install.packages("tmap")
library(tmap)
# Display data
tm_shape(ct) + tm_fill()

# Note the area displayred is too big. 

# Define the area of interest, which can be done by creating a new spatial object, ct_region.
ct_region = st_bbox(c(xmin = -874210.8, xmax = 1051539,
                      ymin = 535177.8, ymax = 2926678),
                    crs = st_crs(ct)) %>% 
  st_as_sfc()

# where did the coordinates for the bounding box come from?  
# we had previously been exploring the "windCF" raster and identified its extent
# get desired spatial extent using bbox()
# load raster library
# library(raster)
# open raster data
# myRas <- raster(x = "week_03/windCF.tif")
# bbox(myRas)
#         min     max
# s1 -874210.8 1051539
# s2  535177.8 2926678

# now that we have created a spatial object representing the extent of the area of interest, we can plot the continent simple feature
ct_map = tm_shape(ct, bbox = ct_region) + tm_fill() 
map2 = ct_map + tm_layout(bg.color = "lightblue")


# Add the Western states
# first create a simple feature object "states"
# use st_read function as before
states <- st_read("week_03/USAboundaries_aea/USAboundaries_aea.shp")
## Reading layer `USAboundaries_aea' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\week_03\USAboundaries_aea\USAboundaries_aea.shp' using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 56 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -4703607 ymin: 456533.4 xmax: 3563081 ymax: 5625473
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
fullWest <- st_read("geogMasks/fullWest.shp")
## Reading layer `fullWest' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\geogMasks\fullWest.shp' using driver `ESRI Shapefile'
## Simple feature collection with 11 features and 57 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -872377.6 ymin: 922080.2 xmax: 1052019 ymax: 2926679
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
map3 = map2 + tm_shape(states) + tm_fill() + 
  tm_shape(fullWest) + tm_polygons(col = "white")
map3

# New shapes can be added with "+ tm_shape(new_obj)". In this case new_obj represents a new spatial object to be plotted on top of preceding layers. 
# The next code chunk uses the function tm_raster() to plot a raster layer (note we can use alpha set to make the layer semi-transparent):

# First we need to import the raster.
# Previously we had used "st_read" to import a shapefile.  
# ct <- st_read("continent/continent.shp")
# This time we will use the raster library to import a raster.
library(raster)
## Loading required package: sp
windSSCat1 <- raster(x = "wind/windCat1CF.tif")
windSSCat2 <- raster(x = "wind/windCat2CF.tif")
windSSCat3 <- raster(x = "wind/windCat3CF.tif")
windSSCat4 <- raster(x = "wind/windCat4CF.tif")

# this is a practice map. 
# map4 = map3 +
#   tm_shape(windSSCat1) + tm_raster(alpha = 0.7)

# Note that the default color palette is yellows and oranges.  That's not right for wind.  Needs to be blue. 
# Use help files to explore color palettes
# help("tm_raster")

# install RcolorBrewer package if you haven't already
# install.packages("RColorBrewer")
# use interactive viewer to choose color palette interactive viewer: http://colorbrewer2.org/
library(RColorBrewer)
#create color ramp variable
ramp <- brewer.pal(6, "Blues")
#decide to use the darker end of the Blues palette
ramp <- c("#9ECAE1", "#6BAED6", "#4292C6", "#2171B5", "#08519C", "#08306B")

map4 = map3 + tm_shape(windSSCat1) + tm_raster(palette = ramp)
map4
## Raster object has 73686898 (9566 by 7703) cells, which is larger than 1e+07, the maximum size determined by the option max.raster. Therefore, the raster will be shown at a decreased resolution of 1e+07 cells. Set tmap_options(max.raster = c(plot = 73686898, view = 73686898)) to show the whole raster.

Now that we have the basic desired elements plotted together in a frame, let’s repeat the frame in a series of facets

solarSSCat1 <- raster(x = "solar/solarCat1CF.tif")
solarSSCat2 <- raster(x = "solar/solarCat2CF.tif")
solarSSCat3 <- raster(x = "solar/solarCat3CF.tif")
solarSSCat4 <- raster(x = "solar/solarCat4CF.tif")
solarCat1Map = map3 + tm_shape(solarSSCat1) + tm_raster(palette = "Oranges") +
   tm_layout(title = "Solar Capacity Factor Category 1") 
solarCat2Map = map3 + tm_shape(solarSSCat2) + tm_raster(palette = "Oranges") +
   tm_layout(title = "Solar Capacity Factor Category 2") 
solarCat3Map = map3 + tm_shape(solarSSCat3) + tm_raster(palette = "Oranges") +
  tm_layout(title = "Solar Capacity Factor Category 3") 
solarCat4Map = map3 + tm_shape(solarSSCat4) + tm_raster(palette = "Oranges") +
  tm_layout(title = "Solar Capacity Factor Category 4") 
tmap_arrange(solarCat1Map, solarCat2Map, solarCat3Map, solarCat4Map)

Repeat the facet chart for wind

windCat1Map = map4 + tm_layout(title = "Wind Capacity Factor Category 1")
windCat2Map = map3 + tm_shape(windSSCat2) + tm_raster(palette = ramp) +
   tm_layout(title = "Wind Capacity Factor Category 2") 
windCat3Map = map3 + tm_shape(windSSCat3) + tm_raster(palette = ramp) +
  tm_layout(title = "Wind Capacity Factor Category 3") 
windCat4Map = map3 + tm_shape(windSSCat4) + tm_raster(palette = ramp) +
  tm_layout(title = "Wind Capacity Factor Category 4") 
tmap_arrange(windCat1Map, windCat2Map, windCat3Map, windCat4Map)

Now we have plotted wind categories 1-4, and solar categories 1-4. The next thing is to create facets for geographies.

fullWest <- st_read("geogMasks/fullWest.shp")
## Reading layer `fullWest' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\geogMasks\fullWest.shp' using driver `ESRI Shapefile'
## Simple feature collection with 11 features and 57 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -872377.6 ymin: 922080.2 xmax: 1052019 ymax: 2926679
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
partWestMask <- st_read("geogMasks/partWestMask.shp")
## Reading layer `partWestMask' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\geogMasks\partWestMask.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 57 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -257148.4 ymin: 1550369 xmax: 1052019 ymax: 2926679
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
inStateMask <- st_read("geogMasks/inStateMask.shp")
## Reading layer `inStateMask' from data source `C:\Users\Emily Leslie\Documents\earth-analytics\data\geogMasks\inStateMask.shp' using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 57 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -856592.1 ymin: 922080.2 xmax: 1052019 ymax: 2926679
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
solarInStateCat1Map = map3 + tm_shape(solarSSCat1) + tm_raster(palette = "Oranges") +
  tm_shape(inStateMask) + tm_fill() + 
  tm_shape(states) +  tm_borders() +
   tm_layout(title = "In-State Solar",legend.position = c("left","bottom"))  

solarPartWestCat1Map = map3 + tm_shape(solarSSCat1) + tm_raster(palette = "Oranges") +
  tm_shape(partWestMask) + tm_fill() + 
  tm_shape(states) +  tm_borders() +
   tm_layout(title = "Part-West Solar",legend.position = c("left","bottom"))  

solarFullWestCat1Map = map3 + tm_shape(solarSSCat1) + tm_raster(palette = "Oranges") +
  tm_shape(states) +  tm_borders() +
   tm_layout(title = "Full-West Solar",legend.position = c("left","bottom")) 

tmap_arrange(solarInStateCat1Map, solarPartWestCat1Map, solarFullWestCat1Map)

windInStateCat1Map = map3 + tm_shape(windSSCat1) + tm_raster(palette = ramp) +
  tm_shape(inStateMask) + tm_fill() + 
  tm_shape(states) +  tm_borders() +
   tm_layout(title = "In-State Wind", legend.position = c("left","bottom")) 

windPartWestCat1Map = map3 + tm_shape(windSSCat1) + tm_raster(palette = ramp) +
  tm_shape(partWestMask) + tm_fill() + 
  tm_shape(states) +  tm_borders() +
   tm_layout(title = "Part-West Wind",legend.position = c("left","bottom"))  

windFullWestCat1Map = map3 + tm_shape(windSSCat1) + tm_raster(palette = ramp) +
  tm_shape(states) +  tm_borders() +
   tm_layout(title = "Full-West Wind",legend.position = c("left","bottom")) 

tmap_arrange(windInStateCat1Map, windPartWestCat1Map, windFullWestCat1Map)

# Use help files to explore layout arguments (fonts, colors, frames, margins)
# help("tm_layout")
# also use help files to explore facet chart options
# help("tm_facets")

next steps: -rearrange so that solar expanding geographies appear in rows, -vertical downward progression toward wider geography -first column: solar capped (instate to fullWest) -second column: wind capped (instate to fullWest) -third column: solar w2w (instate to fullWest progressing down the rows) -fourth column: wind w2w (instate to fullWest progressing down the rows) -legend title should have “capacity factor” spelled out