This RStudio notebook page provides an overview of the code used to plot wind capacity factor maps of the Western US, in the R programming language. This will not address the technical analysis to calculate wind capacity factor from wind speed, nor will it address the analysis to exclude technically infeasible or environmentally sensitive areas from the wind speed data. This page simply focuses on ingesting the spatial data into the R framework, and producing plots with R commands.

When we generate the maps in an R script, this will allow us to create the settings that we want to use, then repeat those same settings for 100 maps, without manually updating the settings for each map, as we would have to do in ArcGIS Pro. With the R script, we can choose the font once for example, and then propagate that same font update through all 100 plots. There are 100 plots because there is one wind energy series for each state, with a subseries for each land use protection level 1-4, and a separate solar subseries for each state, repeated for each land protection level 1-4.

Start by loading libraries. Install the packages if you haven’t already.

# load libraries
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Emily Leslie/Documents/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Emily Leslie/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
# Make sure your working directory is set to  the right file path
# setwd("C:/Users/Emily Leslie/Documents/earth-analytics")

# open raster data
myRas <- raster(x = "data/week_03/windCF.tif")

# plot raster data
plot(myRas,
     main = "West-wide Wind Capacity Factor - Category 1")

# zoom in to one region of the raster
# we may use this later, but need to get the overall plot working first
#plot(myRas,
#  xlim = c(473000, 473030), # define the x limits
#  ylim = c(4434000, 4434030), # define y limits for the plot
#     main = "Lidar Raster - Zoomed into one small region")

Coordinate Reference System

# view coordinate reference system
# assign crs to an object (class) to use for reprojection and other tasks
myCRS <- crs(myRas)
myCRS
## CRS arguments:
##  +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 +ellps=GRS80 +towgs84=0,0,0
#The CRS format, returned by R, is in a PROJ 4 format. This means that the projection information is strung together as a series of text elements, each of which #begins with a + sign.
## +proj=aea The projection of the dataset. Example data are in Albers Equal Area
## +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-114 These are the coordinates of the extent
## +x_0=0
## +y_0=0
# +datum=NAD83 The datum was used to define the center point of the projection. Example raster uses the NAD83 datum.
## +units=m This is the horizontal units that the data are in. Example units are meters.
## +no_defs 
## +ellps=GRS80 
## +towgs84=0,0,0

Extent and Resolution

# what is the x and y resolution for the raster data?
xres(myRas)
## [1] 250
## [1] 250
yres(myRas)
## [1] 250
## [1] 250

# we know from looking at the crs earlier that the linear units are meters
# and I know from working with the data that the cell size is 250 m

What classification values to use?

#show summary of data
summary(myRas)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.14% of all cells)
##             windCF
## Min.    0.04473584
## 1st Qu. 0.30413745
## Median  0.37131137
## 3rd Qu. 0.42378996
## Max.    0.59585685
## NA's    0.00000000
# plot histogram of data
hist(myRas,
     main = "Distribution of raster cell values in the Wind CF data",
     xlab = "Capacity Factor", ylab = "Number of Pixels",
     col = "#6BAED6")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 0% of the raster cells were used. 100000 values used.

Customizing the West-wide wind capacity factor plot (map) Choose color scale, and figure out where to put the breaks in the color scale, based on the quantiles we just identified. For every quantile, place a label on the legend, indicating that capacity factors higher than this value show up as a darker shade of blue on the map.

# 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")

# plot raster data
plot(myRas,
     breaks = c(0.25, 0.30, 0.37, 0.42, 0.59), # note these breaks are identified from the quantiles in the histogram above
     main = "West-wide Wind Capacity Factor - Category 1",
     axes = TRUE,
     col = ramp)

It is nice to see the capacity factorsplotted, but we need state boundaries for reference and orientation. Otherwise these are just pixels floating in space.

Open Vector Layer. The state boundaries are saved as a shape file (vector layer). To open a shapefile you use the readOGR() function. Eventually we are going to use add = TRUE, to add a layer on top of an existing plot in R.

Check the coordinate system for the vector data, in order to get the two layers to line up…

stateBdry <- readOGR("data/week_03/USAboundaries/USAboundaries.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Emily Leslie\Documents\earth-analytics\data\week_03\USAboundaries\USAboundaries.shp", layer: "USAboundaries"
## with 51 features
## It has 54 fields
## Integer64 fields read as strings:  POPULATION POP2010 WHITE BLACK AMERI_ES ASIAN HAWN_PI HISPANIC OTHER MULT_RACE MALES FEMALES AGE_UNDER5 AGE_5_9 AGE_10_14 AGE_15_19 AGE_20_24 AGE_25_34 AGE_35_44 AGE_45_54 AGE_55_64 AGE_65_74 AGE_75_84 AGE_85_UP HOUSEHOLDS HSEHLD_1_M HSEHLD_1_F MARHH_CHD MARHH_NO_C MHH_CHILD FHH_CHILD FAMILIES HSE_UNITS VACANT OWNER_OCC RENTER_OCC
plot(stateBdry,
  main = "Shapefile imported into R - state boundaries",
  axes = TRUE)

crs(stateBdry)
## CRS arguments:
##  +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
## +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

We can see that the raster (wind capacity factors) are in Albers Equal Area projection, and the vector data (state boundaries) are in mercator projection. Use ArcGIS Pro to project the vector data into Abers equal Area projection system, so that the two layers can be overlaid

Once the projections are matched, try plotting again. Remove the axis labels, since they are distracting.

#import the new shapefile in the correct projection
stateBdry <- readOGR("data/week_03/USAboundaries_aea/USAboundaries_aea.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Emily Leslie\Documents\earth-analytics\data\week_03\USAboundaries_aea\USAboundaries_aea.shp", layer: "USAboundaries_aea"
## with 51 features
## It has 56 fields
## Integer64 fields read as strings:  POPULATION POP2010 WHITE BLACK AMERI_ES ASIAN HAWN_PI HISPANIC OTHER MULT_RACE MALES FEMALES AGE_UNDER5 AGE_5_9 AGE_10_14 AGE_15_19 AGE_20_24 AGE_25_34 AGE_35_44 AGE_45_54 AGE_55_64 AGE_65_74 AGE_75_84 AGE_85_UP HOUSEHOLDS HSEHLD_1_M HSEHLD_1_F MARHH_CHD MARHH_NO_C MHH_CHILD FHH_CHILD FAMILIES HSE_UNITS VACANT OWNER_OCC RENTER_OCC
# plot raster data
plot(myRas,
     breaks = c(0.25, 0.30, 0.37, 0.42, 0.59), # note these breaks are identified from the quantiles in the histogram above
     main = "West-wide Wind Capacity Factor - Category 1",
     axes = FALSE,
     col = ramp)
plot(stateBdry,
  main = "Shapefile imported into R - state boundaries",
  add = TRUE)

Looks good. Now all we need to do is refine, and repeat.

Next steps: - add gray continent background (Canada and Mexico)