R Notebook

library(terra)

## terra 1.7.71

library(tidyterra)

## 
## Attaching package: 'tidyterra'

## The following object is masked from 'package:stats':
## 
##     filter

library(sp)
library(ggplot2)
library(ggthemes)
library(urbnmapr)
library(sf)

## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

setwd("C:/Users/konst/Desktop/DIDA 370")


world_map <- rast("NE1_50M_SR_W.tif")
#world_map <- terra::project(world_map, "EPSG:32116", method = "bilinear")


world_map

## class       : SpatRaster 
## dimensions  : 5400, 10800, 3  (nrow, ncol, nlyr)
## resolution  : 0.03333333, 0.03333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : NE1_50M_SR_W.tif 
## colors RGB  : 1, 2, 3 
## names       : NE1_50M_SR_W_1, NE1_50M_SR_W_2, NE1_50M_SR_W_3

plot(world_map)

#map1
setwd("C:/Users/konst/Desktop/DIDA 370/")
statemap <- st_read("states")

## Reading layer `state_map' from data source 
##   `C:\Users\konst\Desktop\DIDA 370\states' using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -125.4117 ymin: 21.73154 xmax: -66.9499 ymax: 49.38436
## Geodetic CRS:  WGS 84

states<-get_urbn_map("states", sf = TRUE)

#filter the data to get just NYS
statemap <- statemap %>% 
  st_transform("EPSG:32116")
print(statemap)

## Simple feature collection with 51 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -5185573 ymin: -1706408 xmax: 1011941 ymax: 2154323
## Projected CRS: NAD83 / New York Central
## First 10 features:
##    state_fips state_abbv  state_name                       geometry
## 1          01         AL     Alabama MULTIPOLYGON (((-853491.4 -...
## 2          04         AZ     Arizona MULTIPOLYGON (((-3440549 -9...
## 3          08         CO    Colorado MULTIPOLYGON (((-2606544 38...
## 4          09         CT Connecticut MULTIPOLYGON (((489885.6 12...
## 5          12         FL     Florida MULTIPOLYGON (((-171491.8 -...
## 6          13         GA     Georgia MULTIPOLYGON (((-574857.7 -...
## 7          16         ID       Idaho MULTIPOLYGON (((-2972226 13...
## 8          18         IN     Indiana MULTIPOLYGON (((-764041.7 -...
## 9          20         KS      Kansas MULTIPOLYGON (((-1935298 32...
## 10         22         LA   Louisiana MULTIPOLYGON (((-940453.2 -...

#ny <- statemap%>% filter(state_name == "NY")

ggplot()+
  geom_spatraster(data = world_map)+
  geom_sf(statemap, mapping = aes(), fill = 'transparent', color = "gray")+
  theme_map()+
  theme_map()+
  scale_fill_whitebox_c(
    palette = "pi_y_g",
    na.value = "white")+
  labs(title = "Landcover in NYS",
       fill = "Fill")+
  theme(plot.title = element_text(hjust = 0.5, size = 16),
        legend.position="right")

## <SpatRaster> resampled to 5e+05 cells for plotting
## ! `tidyterra::geom_spatraster()`: Plotting 3 overlapping layers: NE1_50M_SR_W_1, NE1_50M_SR_W_2, and NE1_50M_SR_W_3. Either:
## 
##   Use `facet_wrap(~lyr)` for faceting or
## 
##   Use `aes(fill = <name_of_layer>)` for displaying single layers

world_map

## class       : SpatRaster 
## dimensions  : 5400, 10800, 3  (nrow, ncol, nlyr)
## resolution  : 0.03333333, 0.03333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : NE1_50M_SR_W.tif 
## colors RGB  : 1, 2, 3 
## names       : NE1_50M_SR_W_1, NE1_50M_SR_W_2, NE1_50M_SR_W_3

usa <- statemap
#%>% filter(state_name == "Alabama")

#The crop function allows us to crop a raster to the extent of a polygon
usa_rast <- crop(world_map, usa)



ggplot()+
  geom_spatraster(data = usa_rast)+
  geom_sf(usa, mapping = aes(), fill = 'transparent', color = "white")+
  theme_map()+
  scale_fill_whitebox_c(
    palette = "pi_y_g",
    na.value = "white")+
  scale_fill_whitebox_c(
    palette = "atlas",
    na.value = "white")+
  labs(title = "Landcover in Broome County",
       fill = "Fill")+
  theme(plot.title = element_text(hjust = 0.5, size = 16),
        legend.position="right")

## <SpatRaster> resampled to 5e+05 cells for plotting
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## ! `tidyterra::geom_spatraster()`: Plotting 3 overlapping layers: NE1_50M_SR_W_1, NE1_50M_SR_W_2, and NE1_50M_SR_W_3. Either:
## 
##   Use `facet_wrap(~lyr)` for faceting or
## 
##   Use `aes(fill = <name_of_layer>)` for displaying single layers

#Install it if you haven't yet
library(leaflet)
library(sf)     
library(dplyr)   

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spData) 

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(urbnmapr)
library(ggplot2)

states2 <- get_urbn_map("counties", sf = TRUE)
states2

## old-style crs object detected; please recreate object with a recent sf::st_crs()

## Warning in CPL_crs_from_input(x): GDAL Message 1: CRS EPSG:2163 is deprecated.
## Its non-deprecated replacement EPSG:9311 will be used instead. To use the
## original CRS, set the OSR_USE_NON_DEPRECATED configuration option to NO.

## old-style crs object detected; please recreate object with a recent sf::st_crs()

#filter the data to get just NYS
nystate <- states2 %>% 
  filter(state_abbv == "NY") %>% 
  st_transform("EPSG:32116")

## old-style crs object detected; please recreate object with a recent sf::st_crs()

plot(nystate)

nystate_data <- read.csv("L-14.csv")
nystate_data

nystate_data$X.1<-NULL
nystate_data$X.2<-NULL
nystate_data$X.4<-NULL
nystate_data$X.5<-NULL
nystate_data$X.7<-NULL
nystate_data$X.8<-NULL
nystate_data$X.10<-NULL
nystate_data$X.11<-NULL
nystate_data$X.13<-NULL
nystate_data$X.14<-NULL
nystate_data$X.16<-NULL
nystate_data$X.17<-NULL
nystate_data$X.19<-NULL
nystate_data$X.21<-NULL
nystate_data$X.22<-NULL
nystate_data$X.23<-NULL
nystate_data$X.24<-NULL
nystate_data$X.25<-NULL

nystate_data2 <- nystate_data[!apply(is.na(nystate_data) | nystate_data == "", 1, all),]

nystate_data2$X.3<- NULL
nystate_data2$X<- NULL
nystate_data2$X.3<- NULL
nystate_data2$X.6<- NULL
nystate_data2$X.9<- NULL
nystate_data2$X.12<- NULL

newdatany<-na.omit(nystate_data2)

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:terra':
## 
##     extract

counties <- nystate %>% 
  separate(county_name, c("county_name", "county"), sep = " County") %>% 
  select(-county)


#GDP..Deflator3

#Now merge it!!
#nycounty<- counties %>%
 # left_join(nystate_data,by = c("county_name" = "county_name"))

#plot(nycounty)

#nycounties.sf = st_transform(nycounty, "EPSG:4326")



p_popup <- paste0("<strong>Inequality: </strong>", newdatany$GDP..Deflator3)

#generate a color palette with five values
#read more about leaflet colors here: https://rstudio.github.io/leaflet/colors.html 
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)

#now plot it with these custom options
leaflet() %>%
  leaflet::addPolygons(data = counties, 
          # remove polygon borders             
                       stroke = FALSE, 
          # set fill color with function from above and value
                       fillColor = ~pal_fun(newdatany$GDP..Deflator3), 
                       fillOpacity = 0.8, 
                       smoothFactor = 0.5, 
          # add pop-up text
                       popup = p_popup)

## Warning: sf layer is not long-lat data

## Warning: sf layer has inconsistent datum (+proj=tmerc +lat_0=40 +lon_0=-76.5833333333333 +k=0.9999375 +x_0=250000 +y_0=0 +datum=NAD83 +units=m +no_defs).
## Need '+proj=longlat +datum=WGS84'