library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(raster)
## Loading required package: sp
library(stars)
## Loading required package: abind
library(tmap)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
source("functions.R")
# converting precipitation rate in kg/m^2/s to total daily inches
pr2inches = function(inpr){
return(inpr * 60*60*24/25.4) }
# converting Kelvin temparature to Fahrenheit
temp2far = function(intemp) {
return(intemp * 9/5 - 459.67) }
# vector version of wgs84 in EPSG
vwgs84 = 4326
# raster version of wgs84 in proj4
rwgs84 = "+proj=longlat +datum=WGS84 +no_defs"
# creating states vector dataframe
usst1 = st_read("cb_2018_us_state_20m.shp") %>%
tolow() %>%
# restricting analysis to coterminous US
filter(name != 'Puerto Rico' &
name != 'American Samoa' &
name != 'United States Virgin Islands' &
name != 'Commonwealth of the Northern Mariana Islands' &
name != 'Alaska' &
name != 'Hawaii' &
name != 'Guam') %>%
st_transform(vwgs84)
## Reading layer `cb_2018_us_state_20m' from data source
## `/Users/apple/Desktop/UNI/FALL 2023/1. CCA/Class/class1923/cb_2018_us_state_20m.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 52 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
## Geodetic CRS: NAD83
# reading infile into stars object and fetching basic metadata
# file is for precipitation from ACCESS-CM2 model,
# SSP3, RCP7.0, years 2075-2100, 1/16th degree resolution,
# monthly data, southeast US, in netcdf format
infile = "pr.ACCESS-CM2.ssp370.r3i1p1f1.2075-2100.LOCA_16thdeg_v20220519.monthly.s_east.nc"
x1 = read_ncdf(infile)
## no 'var' specified, using pr_tavg
## other available variables:
## lon, lat, time
## Will return stars object with 16300440 cells.
## No projection information found in nc file.
## Coordinate variable units found to be degrees,
## assuming WGS84 Lat/Lon.
x1
## stars object with 3 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
## Min. 1st Qu. Median Mean
## pr_tavg [kg/m^2/s] 1.365902e-06 2.62226e-05 4.490687e-05 4.57986e-05
## 3rd Qu. Max. NA's
## pr_tavg [kg/m^2/s] 6.381121e-05 0.0001579557 43651
## dimension(s):
## from to offset delta refsys values
## lon 1 215 -88.75 0.0625 WGS 84 NULL
## lat 1 243 24.69 0.0625 WGS 84 NULL
## time 1 312 NA NA POSIXct 2075-01-16 12:00:00,...,2100-12-16 12:00:00
## x/y
## lon [x]
## lat [y]
## time
# converting to dataframe and remove units (to increase processing speed)
p1 = as_tibble(x1)
units(p1$pr_tavg) = NULL
# converting longitude to traditional +/- longitude
# restricting to rectangle around Georgia
p2 = p1 %>%
mutate(lon = ifelse(lon>180,lon-360,lon)) %>%
filter(lat >= 30 & lat <= 36 &
lon <= -80 & lon >= -86)
# converting daily average precipitation from kg/m^2/s to inches total
# calculating year, month, number of days in month
# calculating precipitation total
p3 = p2 %>%
# filter(year(time) == 2100) %>%
mutate(pr2 = pr2inches(pr_tavg)) %>%
mutate(year = year(time),
month = month(time),
days = days_in_month(time),
pr3 = pr2 * days)
# total precipitation by year
p4 = p3 %>%
group_by(lat,lon,year) %>%
summarize(pr4 = sum(pr3,na.rm=T))
## `summarise()` has grouped output by 'lat', 'lon'. You can override using the
## `.groups` argument.
# calculating average annual precipitation
# removing raster cells with NA
p5 = p4 %>%
group_by(lat,lon) %>%
summarize(pr5 = mean(pr4,na.rm=T)) %>%
select(lon,lat,pr5) %>%
mutate(pr5 = ifelse(pr5==0,NA,pr5))
## `summarise()` has grouped output by 'lat'. You can override using the `.groups`
## argument.
# converting to raster
r5 = rasterFromXYZ(p5,crs=rwgs84)
# convert raster cells to vector polygons
polycells = rasterToPolygons(r5)
# mapping average precipitation for Georgia
tm_basemap("OpenStreetMap") +
tm_shape(r5) +
tm_raster(palette="-viridis") +
tm_shape(usst1) +
tm_borders() +
tm_layout(legend.outside.position = "right",
legend.outside = TRUE) +
tm_shape(polycells) +
tm_borders()
