#################################################
## Processes Temperature data from Perkins paper
## http://www.worldclim.org/current
## Download 2.5 minute average temperature
## Data from 1960-1990
## Place unzipped folder in data/directory (should be named "tmean_2-5m_bil")
## Run script to produce tx_county_temps.csv in data/ folder
## Delete temperature files
#################################################
####################################################
## Processes the historic temperature data
rm(list=ls())
library(tidyverse)
library(rgdal)
library(raster)
extract_county_avg_temp <- function(county, temperature_raster){
county_map <- map_data(map = "county") %>% filter(region=="texas", subregion == county)
if(nrow(county_map)==0){
stop("Problem getting county data, likely spelled a county name incorrectly.")
}
p <- Polygon(coords = county_map[,c("long", "lat")])
ps <- Polygons(list(p), ID = 1)
sps <- SpatialPolygons(list(ps))
mean(extract(temperature_raster, sps)[[1]]/10, na.rm=T)
}
compile_temp_df <- function(raster_location){
## Takes in the raster location and compiles all of the average county temperatures from that raster
temperature_raster <- raster(raster_location) # "data/tmean_2-5m_bil/tmean1.bil")
counties <- map_data(map="county") %>% filter(region=="texas") %>% dplyr::select(subregion) %>% unique()
county_temps <- purrr::map(counties$subregion, extract_county_avg_temp, temperature_raster)
counties$avg_temp <- unlist(county_temps)
## Finds the number of the file being called, and sets to that month abbreviation
counties$month <- month.abb[as.numeric(strsplit(strsplit(raster_location, split = ".", fixed = T)[[1]], split="n", fixed=T)[[1]][3])]
counties
}
raster_locations <- list.files(path="data/tmean_2-5m_bil", pattern = "*.bil", full.names = T)
tx_temperatures <- purrr::map(raster_locations, compile_temp_df)
tx_temperatures <- tx_temperatures %>% bind_rows() %>% mutate(month = factor(month, levels = month.abb)) %>%
arrange(month)
write_csv(tx_temperatures, path = "data/tx_historic_temps.csv")
#########################################################
## Processes the actual temperature data from 2016-17
## 2017 data: https://www.esrl.noaa.gov/psd/repository/entry/show?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3AL05BUlIvRGFpbGllcy9tb25vbGV2ZWwvYWlyLjJtLjIwMTcubmM%3D
## 2016 data: https://www.esrl.noaa.gov/psd/repository/entry/show?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3AL05BUlIvRGFpbGllcy9tb25vbGV2ZWwvYWlyLjJtLjIwMTYubmM%3D
library(ncdf4)
library(lubridate)
t_2016 <- brick("data/air.2m.2016.nc")
t_2017 <- brick("data/air.2m.2017.nc")
get_county_temps <- function(current_county, brick_data, current_year){
## Extract the specified county spatial information
county_map <- map_data(map = "county") %>% filter(region=="texas", subregion == current_county)
## First turn the county into a spatial object on the correct projection
p <- Polygon(coords = county_map[,c("long", "lat")])
ps <- Polygons(list(p), ID = 1)
sps <- SpatialPolygons(list(ps), proj4string = CRS("+proj=longlat"))
sps <- spTransform(sps, projection(t_2016))
## Then extract the temperature data
cty_temps <- extract(brick_data, sps)
## If list is longer than 1 need to make sure rest properly works
if(length(cty_temps)!=1){
browser()
}
## Average all readings within county for each date
cty_temps <- apply(cty_temps[[1]], MARGIN = 2, FUN = mean)
## Convert to data_frame with proper units, and extract and average across month
data_frame(date = names(cty_temps), temperature = cty_temps-273.15) %>%
separate(col = date, sep = "\\.", into = c("year", "month", "day", "hour", "minute", "second")) %>%
mutate(year = substring(year, 2)) %>%
unite(col = date, year:day, sep="-") %>%
mutate(date = ymd(date)) %>%
filter(year(date) == current_year) %>%
group_by(month = month(date, label=TRUE)) %>%
summarise(avg_temp = mean(temperature)) %>%
mutate(subregion = current_county, year = current_year) %>%
dplyr::select(subregion, year, month, avg_temp)
}
## extract counties
counties <- map_data(map = "county") %>% filter(region=="texas") %>%
group_by(subregion) %>%
summarize(county=unique(subregion)) %>%
ungroup() %>%
dplyr::select(county)
## Run for all counties of 2016
county_temps_2016 <- counties[[1]] %>% purrr::map(~get_county_temps(current_county =.x, brick_data=t_2016, current_year=2016)) %>% bind_rows()
county_temps_2017 <- counties[[1]] %>% purrr::map(~get_county_temps(current_county =.x, brick_data=t_2017, current_year=2017)) %>% bind_rows()
write_csv(bind_rows(county_temps_2016, county_temps_2017), path = "data/tx_actual_temps.csv")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.