class: Environmental Data Analysis title: "Assignment 2 Lab 2" author: "Riley Chan (rhc2157)" output: html_document ---

Using same project as Assignment 1- so data folder should be similar

load libraries

library(tidyverse) library(readr) library(janitor) library(purrr) library(repurrrsive) library(ggplot2) library(htmltools) library(dplyr) library(sf) library(wesanderson) library(sf) library(leaflet) library(gstat) library(stars) library(tmap) ***

QUESTION 1 WORK

```

Get the list of available palettes

wespalettes <- wespalettes

Calculate palette lengths

palettelength <- mapdbl(wes_palettes, ~length(.))

Filter palettes with more than 4 colors

filterpalette <- wespalettes[palette_length > 4]

Get the names of the selected palettes

wesSelectPals <- c(filter_palette) head(wesSelectPals) typeof(wesSelectPals) ```

QUESTION 2 WORK

```

Load cs. data from Assignment 1

imt<- read_csv("Data1/Trade/InternationalMerchandiseTrade.csv", skip = 1)

Clean and filter data to get Balance imports/exports

imt<- clean_names(imt) imt<- imt%>% select(x2, year, series, value) %>% filter(series=='Balance imports/exports (millions of US dollars)')%>% rename(Country = x2)

Get average from sorted data

imtaverage <- imt[205:1680,] %>% groupby(Country)%>% summarize(averagebalance = mean(value)) %>% arrange(desc(averagebalance))%>% top_n(5)

Arranges data for top 5 countries

topfive <- imt[205:1680,] %>% filter(Country %in% imtaverage$Country, year %in% year)%>% group_by(Country)

Create a ggplot line graph using the current palette

Got stuck here for a while- might not be the right answer

Used for statement to go through all color palettes in wesSelectPals

Cannot figure out how to get print the name of each palette in the title

Use for loop to iterate color palettes in list of wesSelectPals

for (i in wesSelectPals){ mapplot<- ggplot(topfive, aes(x = year, y = value, color= Country)) + geomline() + scalecolormanual(values = i)+ labs( title = paste("Average of Balance of Imports/Exports in",wesSelectPals, "Palette"), x = "Year", y = "Balance Imports/Exports", name = "Country" ) print(mapplot) } ```

QUESTION 3 WORK

```

Load in cb_2018 shapefile and ACS data

cb2018 <- stread("Data2/cb2018uscounty20m/cb2018uscounty20m.shp") acs_2018 <- read.csv("Data2/ACS2018Counties.csv")

Rename Geo_FIPS column in GEOID in order to join datasets

acs2018<- acs2018%>% rename(GEOID= FIPS)

Transform cb_2018 shapefile to north american albers equal area projection

cb2018$geometry <-sttransform(cb2018$geometry, "+proj=aea +lat1=20 +lat2=60 +lat0=40 +lon0=-96 +x0=0 +y0=0 +ellps=GRS80 +datum=NAD83 +units=mi +nodefs")

After renaming Geo_FIPS, merge shapefile and ACS data

cbacs2018merge <- merge(cb2018, acs2018, by = "GEOID", all = TRUE ) head(cbacs2018merge)

Calculate areaCalc of geometry in square miles

areaCalc <- starea(cbacs2018merge$geometry)

Adds areaCalc as column

cbacs2018_merge$areaCalc<- areaCalc

Print areaCalc

print(areaCalc)

Cleans headers

cbacs2018merge<- cleannames(cbacs2018_merge)

Creates new column for population density using total population and the calculated area

cbacs2018merge$popDensity <- as.integer(cbacs2018merge$total_population)/ areaCalc

Print population density

print(cbacs2018_merge$popDensity) ```

QUESTION 4 WORK

```

I cannot figure this out I keep getting different error messages. At wits end on this question

Been stuck on it for a few hours and cannot get this to work

Keep getting this error "Error in to_ring.default(x) :

Don't know how to get polygon data from object of class XY,GEOMETRYCOLLECTION,sfg"

cbacs2018merge <- sttransform(stassf(cbacs2018_merge), 4326)

quantilesmerge <- colorQuantile(palette = "YlOrRd", domain = cbacs2018merge$popDensity, n = 5)

leafletmap <- leaflet(cbacs2018merge)%>% addTiles()%>% addPolygons(stroke = FALSE, color = ~quantiles_merge(popDensity), fillOpacity = .7, label = ~paste(NAME, "
", format(popDensity, digits =2)), labelOptions= labelOptions(direction = "auto")%>% setView(lat=43.61, lng= -116.5, zoom = 3))

Error seems to be with the addPolygons line- NO IDEA WHAT TO DO! Spent a few hours trying to research this, had no luch

print(leaflet_map) ```

Question 5 WORK

```

Loads air quality .csv

airquality <- read.csv("Data2/advizplotvaldata.csv")

Creates new column for statefp with code "06"

air_quality$statefp <- "06"

Clean data headers

airquality<- cleannames(air_quality)

maxpm25data<- airquality%>% groupby(county, statefp) %>% summarize(maxpm25 = max(dailymeanpm25_concentration))

centroidmerge <- innerjoin(cbacs2018merge, maxpm25data, by = c('name' = 'county')) print(centroid_merge)

centroidmergeca<- centroidmerge%>% filter(statefips=='06') print(centroidmergeca)

centroidmergeca$centroid <- stcentroid(centroidmerge_ca)

Comment inner join removes values that are not the same

Prior to inner join, I merged the data using left_join

However, left_join gave different values, so I used inner join insead

```

QUESTION 6 WORK

```

Cast the centroid merged data as Spatial

centroidmergecasp <- asSpatial(centroidmergeca, cast = TRUE)

sfgrid<- stmakegrid(centroidmerge_ca, n= 150)

centroidsfcaidw <- idw(centroidmergeca$maxpm25 ~1, centroidmergecasp, sf_grid, idp=2)

Rasterized the idw data using st_rasterize with the field var1.pred

aqraster <- strasterize(centroidsfcaidw, field = "var1.pred")

Map the raster data using tmap

tmapraster <- tmshape(aqraster) + tmraster("var1.pred", style = "quantile", n = 9, palette = colorRampPalette(c("lightgoldenrod", "firebrick3"))(n = 9))+ tmlayout(bg.color = "white", legend.outside = TRUE)+ tmshape(centroidmergeca) + tmborders(col = "black", lwd = 1, lty = "solid")+ tmshape(centroidmergeca$centroid)+ tmsymbols(col = "grey8", size = .15)

Print map

print(tmap_raster) ```

QUESTION 7 WORK

```

Convert rasterdata to sf

aqrastersf <- stas_sf(aqraster)

Creates new field with st_intersecton to find the ovrlapping values

aqrasterintersect <- stintersection( aqrastersf, centroidmergeca)

Creates new data intersection with aggregate by list column "geoid" with function max.

maxpm25intersect <- aggregate(aqrasterintesect["var1.pred"], by = list(aqraster_intersect$geoid), FUN = max)

Create new field to give values of var1.pred that are over 200

over200 <- maxpm25intersect[maxpm25_intersect$var1.pred >200,]

over_200

Getting the Error "Error: object 'aqrasterintersect' not found"- might be pulling from the wrong data set.

My computer also became really slow when loading the code at the end of the document- need to figure out how to keep usability as code gets longer

```