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
113
114
115
|
R version 3.6.2 (2019-12-12) -- "Dark and Stormy Night"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> ###################################################
> ### chunk number 19:
> ###################################################
> suppressPackageStartupMessages(library(sp))
> suppressPackageStartupMessages(library(spacetime))
> if (require(gstat)) {
+ data(wind)
+ wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]])))
+ wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]])))
+ coordinates(wind.loc) = ~x+y
+ proj4string(wind.loc) = "+proj=longlat +datum=WGS84"
+
+
+ ###################################################
+ ### chunk number 20:
+ ###################################################
+ suppressPackageStartupMessages(library(mapdata))
+ plot(wind.loc, xlim = c(-11,-5.4), ylim = c(51,55.5), axes=T, col="red")
+ map("worldHires", add=T, col = grey(.5))
+ text(coordinates(wind.loc), pos=1, label=wind.loc$Station, cex=.7)
+
+
+ ###################################################
+ ### chunk number 21:
+ ###################################################
+ wind[1:3,]
+ wind.loc[1:3,]
+
+
+ ###################################################
+ ### chunk number 22:
+ ###################################################
+ wind$time = ISOdate(wind$year+1900, wind$month, wind$day)
+ wind$jday = as.numeric(format(wind$time, '%j'))
+
+
+ ###################################################
+ ### chunk number 23:
+ ###################################################
+ # match order of columns in wind to Code in wind.loc;
+ # convert to utm zone 29, to be able to do interpolation in
+ # proper Euclidian (projected) space:
+ pts = coordinates(wind.loc[match(names(wind[4:15]), wind.loc$Code),])
+ pts = SpatialPoints(pts)
+ proj4string(pts) = "+proj=longlat +datum=WGS84"
+ suppressPackageStartupMessages(library(rgdal))
+ pts = spTransform(pts, CRS("+proj=utm +zone=29 +datum=WGS84"))
+ stations = 4:15
+ # note the t() in:
+ w = STFDF(pts, wind$time, data.frame(values = as.vector(t(wind[stations]))))
+
+ suppressPackageStartupMessages(library(maptools))
+ m = map2SpatialLines(
+ map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5), plot=F))
+ proj4string(m) = "+proj=longlat +datum=WGS84"
+ m = spTransform(m, CRS("+proj=utm +zone=29 +datum=WGS84"))
+
+ # setup grid
+ grd = SpatialPixels(SpatialPoints(makegrid(m, n = 300)),
+ proj4string = proj4string(m))
+ # grd$t = rep(1, nrow(grd))
+ #coordinates(grd) = ~x1+x2
+ #gridded(grd)=TRUE
+
+ # select april 1961:
+ w = w[, "1961-04"]
+
+ covfn = function(x,y) {
+ du = spDists(coordinates(x), coordinates(y))
+ t1 = as.numeric(index(x)) # time in seconds
+ t2 = as.numeric(index(y)) # time in seconds
+ dt = abs(outer(t1, t2, "-"))
+ # separable, product covariance model:
+ 0.6 * exp(-du/750000) * exp(-dt / (1.5 * 3600 * 24))
+ }
+
+ n = 10
+ tgrd = seq(min(index(w)), max(index(w)), length=n)
+ if (FALSE) {
+ pred = krige0(sqrt(values)~1, w, STF(grd, tgrd), covfn)
+ wind.pr = STFDF(grd, tgrd, data.frame(pred = pred))
+
+
+ ###################################################
+ ### chunk number 24:
+ ###################################################
+ spl = list(list("sp.points", pts, first=F, cex=.5),
+ list("sp.lines", m, col='grey'))
+ stplot(wind.pr, col.regions=bpy.colors(),
+ par.strip.text = list(cex=.5), sp.layout = spl)
+ summary(wind.pr)
+ }
+ }
Loading required package: gstat
>
> proc.time()
user system elapsed
1.061 0.073 1.127
|
