-
Notifications
You must be signed in to change notification settings - Fork 0
/
interpolation.R
67 lines (49 loc) · 2.09 KB
/
interpolation.R
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
###########################################
#Interpolation
library(lubridate)
library(tidyverse)
library(vegan)
library(sf)
if (!require("rspatial")) devtools::install_github('rspatial/rspatial')
library(rspatial)
library( spgwr )
library(gstat)
#load the data
temps = readRDS("Temp_filtered (1).rds")
#first calculate the daily means, mins, and maxes
tempminmax = temps %>%
filter(Date > as.Date("2000-1-1"), Station != "RYF") %>%
mutate(julian = yday(Date)) %>%
group_by(Station, julian) %>%
summarize(Temp = mean(Temp, na.rm = T), Tempmax = max(Temp, na.rm = T),
Tempmin = min(Temp, na.rm = T))
#read in shapefile of the delta
delta = read_sf("deltashap/hydro_delta_marsh.shp")
#add lat/longs for the stations
stas = read.csv("StationLatLongs.csv")
#attached lat/longs to mean temperature
temp2 = merge(tempminmax, stas)
#Specify a coordinate reference system and turn it into a spatial object
alb <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
#Filter the Delta shapefile so it doesn't include areas where we dont have a lot of data
delta2 = sf::as_Spatial(dplyr::filter(delta, HNAME != "SAN FRANCISCO BAY", HNAME != "SAN PABLO BAY"))
#I'm not positive we have to do this transformation, but it was in teh example so I'm doing it.
ctst <- spTransform(delta2, alb)
# just look at summer temps
temp2.1 = dplyr::filter(temp2, julian < 280, julian >180)
sp = temp2.1
coordinates(sp) = ~ Longitude + Latitude
crs(sp) <- "+proj=longlat +datum=NAD83"
spt <- spTransform(sp, alb)
#Raster of the Delta
r <- raster(ctst, res=100)
vr <- rasterize(vca, r, 'prec')
gs <- gstat(formula=Tempmax~1, locations=spt)
idw <- interpolate(r, gs)
## [inverse distance weighted interpolation]
idwr <- mask(idw, ctst)
plot(idwr, xlim = c(-200000, -100000), ylim = c(-50000,50000))
###############################################
#look at dat from those weird stations
test = filter(temps, Station %in% c("GSS", "NMR", "SWE", "GES", "SMR"))
ggplot(test, aes(x = Datetime, y = Temp)) + geom_line(aes(color = Station)) + facet_wrap(~Station)