load("wind.RData") require(sp) summary(wind) wind.loc 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 # fig 1: if (require(mapdata)) { map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5)) plot(wind.loc, add=TRUE, pch=16) text(coordinates(wind.loc), pos=1, label=wind.loc$Station) } wind$time = ISOdate(wind$year+1900, wind$month, wind$day) # time series of e.g. Dublin data: plot(DUB~time, wind, type= 'l', ylab = "windspeed (knots)", main = "Dublin") # fig 2: #wind = wind[!(wind$month == 2 & wind$day == 29),] wind$jday = as.numeric(format(wind$time, '%j')) windsqrt = sqrt(0.5148 * wind[4:15]) Jday = 1:366 daymeans = apply(sapply(split(windsqrt - mean(windsqrt), wind$jday), mean), 2, mean) plot(daymeans ~ Jday) lines(lowess(daymeans ~ Jday, f = 0.1)) # subtract the trend: meanwind = lowess(daymeans ~ Jday, f = 0.1)$y[wind$jday] velocity = apply(windsqrt, 2, function(x) { x - meanwind }) # match order of columns in wind to Code in wind.loc: pts = coordinates(wind.loc[match(names(wind[4:15]), wind.loc$Code),]) # fig 3, but not really yet... dists = apply(pts, 1, function(x) spDistsN1(pts, x, longlat=TRUE)) corv = cor(velocity) sel = !(as.vector(dists) == 0) plot(as.vector(corv[sel]) ~ as.vector(dists[sel]), xlim = c(0,500), ylim = c(.4, 1), xlab = "distance (km.)", ylab = "correlation") # plots all points twice, and zero distance # now really get fig 3: ros = rownames(corv) == "ROS" dists.nr = dists[!ros,!ros] corv.nr = corv[!ros,!ros] sel = !(as.vector(dists.nr) == 0) plot(as.vector(corv.nr[sel]) ~ as.vector(dists.nr[sel]), pch = 3, xlim = c(0,500), ylim = c(.4, 1), xlab = "distance (km.)", ylab = "correlation") # add outlier: points(corv[ros,!ros] ~ dists[ros,!ros], pch=16, cex=.5) xdiscr = 1:500 lines(xdiscr, .968 * exp(- .00134 * xdiscr))