library(WikipediR)
library(tidyverse)

content=page_content(
  page_name = "API",
  project="Wikipedia",
  language="en",
  as_wikitext=T)

c1=content$parse$wikitext%>%
  str_split(boundary("sentence"),simplify = T)%>%
  str_replace_all("'","")%>%
  str_replace_all("\\[|\\]","")
# results:
cat(c1[5:6])

library(tidyverse)
library(httr)
library(jsonlite)

library(tidycensus)

net_migration <- get_estimates(geography = "county",
                               variables = "RNETMIG",
                               geometry = TRUE) %>% 
                               tigris::shift_geometry()


net_migration <- net_migration %>%
  mutate(groups = case_when(
    value > 15 ~ "+15 and up",
    value > 5 ~ "+5 to +15",
    value > -5 ~ "-5 to +5",
    value > -15 ~ "-15 to -5",
    TRUE ~ "-15 and below"
  ))

ggplot() +
  geom_sf(data = net_migration, aes(fill = groups, color = groups), lwd = 0.1) +
  geom_sf(data = tidycensus::state_laea, fill = NA, color = "black", lwd = 0.1) +
  scale_fill_brewer(palette = "PuOr", direction = -1) +
  scale_color_brewer(palette = "PuOr", direction = -1, guide = "none") +
  coord_sf() +
  labs(title = "Net migration per 1000 residents by county",
       subtitle = "US Census Bureau 2017 Population Estimates",
       fill = "Rate",
       caption = "Data acquired with the R tidycensus package")

library(rnassqs)
#nassqs_auth("<your api key>")

df <- nassqs_acres(
  commodity_desc = "CORN",
  year=2000:2024,
  county_name = "Erie",
  agg_level_desc = "COUNTY",
  state_alpha = c("NY"),
  statisticcat_desc = "AREA") %>% 
  filter(short_desc=="CORN - ACRES PLANTED") %>% 
  dplyr::select(year, Value) %>%
  transmute(
    year = as.numeric(gsub(",", "", year)),
    value = as.numeric(gsub(",", "", Value)))

geocode("University at Buffalo, NY")

geocode("Washington")

geocode("Washington D.C.")

UB=geocode("University at Buffalo, Buffalo, NY") |>
  st_as_sf(coords=c(1,2)) |> #convert to sf
  st_set_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") |>
  st_transform(32617) |> #transform to UTM
  st_buffer(dist= 20000) |>  #radius of buffer in decimal degrees
  st_transform(4326) #transform back to longlat

library(leaflet)
leaflet() %>%
  addTiles() %>%
  addPolygons(data=UB) %>% 
  addMarkers(data=st_centroid(UB))

# devtools::install_github("ropensci/FedData")
library(FedData)
library(terra)
library(tidyverse)

# Get the NED (USA ONLY)
# Returns a raster
NED <- get_ned(
  template = UB,
  label = "UB",res = 13
)

library(rasterVis)

gplot(NED)+
  geom_raster(aes(fill=value))+
  scale_fill_viridis_c()+
  geom_sf(data=UB,inherit.aes = F,fill="transparent",col="red")+
  geom_sf(data=st_centroid(UB),inherit.aes = F,fill="transparent",col="red")

# Get the DAYMET (North America only)
# Returns a raster
DAYMET <- get_daymet(
  template = UB,
  label = "UB",
  elements = c("tmax"),
  years = 2023:2023
) 

gplot(DAYMET$tmax[[1:10]])+
  geom_raster(aes(fill=value))+
  facet_wrap(~variable)+
  scale_fill_viridis_c()+
  geom_sf(data=st_transform(UB,crs(DAYMET$tmax)),inherit.aes = F,fill="transparent",col="red")+
  labs(title="Daymet Maximum Temperature")

# Get the daily GHCN data (GLOBAL)
# Returns a list: the first element is the spatial locations of stations,
# and the second is a list of the stations and their daily data
GHCN.tmax <- get_ghcn_daily(
  template = UB,
  label = "UB",
  elements = c("tmax")
)

tmax <- GHCN.tmax$spatial %>%
  dplyr::mutate(label = paste0(ID, ": ", NAME))

leaflet() %>%
  addTiles() %>%
  addMarkers(data=tmax)

kable(head(GHCN.tmax$tabular$USW00014733$TMAX))

daily <- GHCN.tmax$tabular$USW00014733$TMAX |>
  gather(day,value,-STATION,-YEAR,-MONTH) |>
  mutate(day=str_remove(day,"D"),
         tmax=value/10,
         date=as_date(paste(YEAR,MONTH,day,sep="-"))) %>% 
  arrange(date)

kable(head(daily))

ggplot(daily,
       aes(y=tmax,x=date))+
  geom_line(col="red")

daily_recent=filter(daily,date>as.Date("2020-01-01"))

ggplot(daily_recent,
        aes(y=tmax,x=date))+
 geom_line(col="red")

library(zoo)
d_rollmean = daily_recent %>% 
  arrange(date) %>%
  mutate(tmax.60 = rollmean(x = tmax, 60, align = "center", fill = NA),
         tmax.b60 = rollmean(x = tmax, 60, align = "right", fill = NA))

d_rollmean%>%
  ggplot(aes(y=tmax,x=date))+
    geom_line(col="grey")+
    geom_line(aes(y=tmax.60),col="red")+
    geom_line(aes(y=tmax.b60),col="darkred")

ggplot(daily_recent,aes(y=tmax,x=lag(tmax)))+
  geom_point()+
  geom_abline(intercept=0, slope=1)

acf(tmax.ts,lag.max = 365*3,na.action = na.exclude )

pacf(tmax.ts,lag.max = 365,na.action = na.exclude )

1938

round(1938,-1)

floor(1938/10)

floor(1938/10)*10 

daily2=daily%>%
  mutate(month=as.numeric(format(date,"%m")),
        year=as.numeric(format(date,"%Y")),
        season=case_when(
          month%in%c(12,1,2)~"Winter",
          month%in%c(3,4,5)~"Spring",
          month%in%c(6,7,8)~"Summer",
          month%in%c(9,10,11)~"Fall"),
        dec=(floor(as.numeric(format(date,"%Y"))/10)*10))
daily2%>%dplyr::select(date,month,year,season,dec,tmax)%>%head()%>%kable()

daily2%>%
  mutate(period=ifelse(year<=1976-01-01,"early","late"))%>% #create two time periods before and after 1976
  group_by(period)%>%  # divide the data into the two groups
  summarize(n=n(),    # calculate the means between the two periods
            tmax=mean(tmax,na.rm=T))%>%
  kable()

daily2%>%
  group_by(year)%>%
  summarize(n=n())%>%
  ggplot(aes(x=year,y=n))+
  geom_line()

# which years don't have complete data?
daily2%>%
  group_by(year)%>%
  summarize(n=n())%>%
  filter(n<360)

daily2%>%
  filter(year>1938, year<2024)%>%
  group_by(dec)%>%
  summarize(
            n=n(),
            tmax=mean(tmax,na.rm=T)
            )%>%
  ggplot(aes(x=dec,y=tmax))+
  geom_line(col="grey")

df=daily2%>%
  mutate(doy=as.numeric(format(date,"%j")),
         doydate=as.Date(paste("2024-",doy),format="%Y-%j"))

ggplot(df,aes(x=doydate,y=tmax,group=year))+
  geom_line(col="grey",alpha=.5)+ # plot each year in grey
  stat_smooth(aes(group=1),col="black")+   # Add a smooth GAM to estimate the long-term mean
  geom_line(data=filter(df,year==2024),col="red")+  # add 2023 in red
  scale_x_date(labels = date_format("%b"),date_breaks = "2 months") + 
  xlab("Day of year")

ggplot(df,aes(x=doydate,y=tmax,group=year))+
  geom_line(col="grey",alpha=.5)+
  stat_smooth(aes(group=1),col="black")+
  geom_line(data=filter(df,year==2024),col="red")+
  scale_x_date(labels = date_format("%b"),date_breaks = "2 months",
               lim=c(as.Date("2024-08-01"),as.Date("2024-10-31")))

seasonal=daily2%>%
  group_by(year,season)%>%
  summarize(n=n(),
            tmax=mean(tmax))%>%
  filter(n>75)

seasonal |>
    filter(!is.na(season)) |>
ggplot(aes(y=tmax,x=year))+
  facet_wrap(~season,scales = "free_y")+
  stat_smooth(method="lm", se=T)+
  geom_line()

lm1=seasonal%>%
  filter(season=="Summer")%>%
  lm(tmax~year,data=.)

summary(lm1)$r.squared

tidy(lm1)%>%kable() 

# fit a lm model for each group
  daily2 %>%
  filter(!is.na(season)) |>
  group_by(season) %>% #process by season
  nest() %>% # cut into groups
  mutate(
    lm_tmax = purrr::map(data, .f = ~ lm(tmax ~ year, data = .)), #fit model for each group
    tmax_tidy = purrr::map(lm_tmax, broom::tidy) #summarize model for each group
  ) %>%
  unnest(tmax_tidy) %>%
  filter(term == "year") %>% 
  dplyr::select(-data, -lm_tmax) %>% 
  kable()