WindClimate/Oneyeartest.R at main · rmgriffin/WindClimate · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
## Use this to start every program.  This clears out previous information from memory
rm(list=ls())

## Initalize renv for library lockfile
library(renv)
#renv::init()

## Packages
#Sys.setenv(RENV_PATHS_RTOOLS = "C:/rtools40/") # https://github.com/rstudio/renv/issues/225

PKG <- c("raster","sf","tidyverse","rgdal","ncdf4","RColorBrewer","lattice","googledrive","tmap","chron","furrr","nngeo")

for (p in PKG) {
  if(!require(p,character.only = TRUE)) {
    install.packages(p)
    require(p,character.only = TRUE)}
}
rm(p,PKG)

## Snapshot of libraries used
renv::snapshot()

## Downloading supporting data
# Download from google drive to directory "Data"
setwd("~/Github/WindClimate")
dir.create(file.path('Data'), recursive = TRUE)
folder_url<-"https://drive.google.com/open?id=1Yc9RiqXXj0qja2DaQaXBwyi267mMZ4Yd"
folder<-drive_get(as_id(folder_url))
files<-drive_ls(folder)
dl<-function(files){
  walk(files, ~ drive_download(as_id(.x), overwrite = TRUE))
}
setwd("./Data")
system.time(map(files$id,dl))
setwd("..")
rm(files, folder, folder_url, dl)

## Download file https://na-cordex.org/data.html
options(timeout=100000) # Keeps download.file from returning a timeout for long-downloading files
download.file("ftp://anonymous:anonymous@ftp2.psl.noaa.gov/pub/Public/dswales/wrfout_d01_2007.nc", destfile = "./Data/wrfout_d01_2007.nc",)

## Actions with a NetCDF
# Reading data
df<-nc_open("./Data/wrfout_d01_2007.nc")
#df1<-nc_open("./Data/sfcWind.hist.GFDL-ESM2M.RegCM4.day.NAM-22i.raw.nc")
# dfb<-brick("./Data/wrfout_d01_2007.nc", crs = 4326)
# dfb<-projectRaster(dfb, crs = 3857)
cst<-st_read("./Data/global_polygon.gpkg")
cst<-st_transform(cst,crs = st_crs(3857))
eeza<-st_read("./Data/eez_atlantic.gpkg") # US EEZ vector https://www.marineregions.org/gazetteer.php?p=details&id=8456
eeza<-st_transform(eeza,st_crs(3857))
eeza<-eeza[1:2]
print(df)

# p1<-tm_shape(dfb[[1003]]) +
#   tm_raster(palette = "Greens", colorNA = NULL, title = "Wind Speed (m/s)") +
#   # tm_shape(cst) +
#   # tm_borders(col = "black", lty = "dashed") +
#   tm_layout(legend.outside = TRUE)
# p1

# Exploring variables
y<-ncvar_get(df,"lat")
ny<-dim(y)
head(y)

x<-ncvar_get(df,"lon")
nx<-dim(x)
head(x)

print(c(nx,ny))

year<-ncvar_get(df,"year")
head(year)

day<-ncvar_get(df,"day")
head(day)

hour<-ncvar_get(df,"hour") # One observation every 3 hours
head(hour)

z<-ncvar_get(df,"level")
head(z)

# Getting data
wsu<-as.vector(ncvar_get(df,"u"))
head(wsu)
wsv<-as.vector(ncvar_get(df,"v"))
head(wsv)

# Domain
x2<-as.matrix(expand.grid(x))
y2<-as.matrix(expand.grid(y))
xy<-as.data.frame(cbind(x2,y2))
names(xy)<-c("lon","lat")
xy<-st_as_sf(xy, coords = c("lon", "lat"),crs=4269, remove = FALSE)
xy<-st_transform(xy,st_crs(3857))

ggplot() + # Plot of observation locations within domain
  geom_sf(data = cst) +
  geom_sf(data = xy)

# Creating a dataframe from whole array
ws<-sqrt(wsu^2+wsv^2)
hist(ws, xlim = c(0, 100), breaks = seq(0,max(ws)+2,2))
max(ws)
length(ws[ws > 32.7])/length(ws) # Ratio of observations greater than hurricane speed
wsday<-colMeans(matrix(ws, nrow=8)) # A vector of daily mean wind speeds, summarize from a vector of 3 hr wind speeds
expand.grid.df<-function(...) Reduce(function(...) merge(..., by=NULL), list(...)) # https://stackoverflow.com/questions/11693599/alternative-to-expand-grid-for-data-frames
system.time(xyz<-expand.grid.df(xy,z,seq(1,365,1))) # Creating a long format dataframe that is # of grid points * # of hub heights * # of days
xyz$wsday<-wsday
head(xyz)