Custom Query (1118 matches)

Action item for release

We need to update the presentation with "what is new in 12.0" and make sure the correct applications show up in the slides.

due to a missing dependency maps

> library(gstat)
> demo(wind)


	demo(wind)
	---- ~~~~

Type  <Return>	 to start : 

> #pdf("wind.pdf")
> library(gstat)

> library(rgdal)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.1, released 2012/05/15
Path to GDAL shared files: /usr/share/gdal/1.9
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files: (autodetected)

> library(maptools)
Loading required package: foreign
Checking rgeos availability: FALSE
 	Note: when rgeos is not available, polygon geometry 	computations in maptools depend on gpclib,
 	which has a restricted licence. It is disabled by default;
 	to enable gpclib, type gpclibPermit()

> # load wind data, run test:
> example(wind)

wind> data(wind)

wind> summary(wind)
      year          month             day             RPT       
 Min.   :61.0   Min.   : 1.000   Min.   : 1.00   Min.   : 0.67  
 1st Qu.:65.0   1st Qu.: 4.000   1st Qu.: 8.00   1st Qu.: 8.12  
 Median :69.5   Median : 7.000   Median :16.00   Median :11.71  
 Mean   :69.5   Mean   : 6.523   Mean   :15.73   Mean   :12.36  
 3rd Qu.:74.0   3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:15.92  
 Max.   :78.0   Max.   :12.000   Max.   :31.00   Max.   :35.80  
      VAL             ROS             KIL              SHA       
 Min.   : 0.21   Min.   : 1.50   Min.   : 0.000   Min.   : 0.13  
 1st Qu.: 6.67   1st Qu.: 8.00   1st Qu.: 3.580   1st Qu.: 6.75  
 Median :10.17   Median :10.92   Median : 5.750   Median : 9.96  
 Mean   :10.65   Mean   :11.66   Mean   : 6.306   Mean   :10.46  
 3rd Qu.:14.04   3rd Qu.:14.67   3rd Qu.: 8.420   3rd Qu.:13.54  
 Max.   :33.37   Max.   :33.84   Max.   :28.460   Max.   :37.54  
      BIR              DUB              CLA              MUL        
 Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
 1st Qu.: 4.000   1st Qu.: 6.000   1st Qu.: 5.090   1st Qu.: 5.370  
 Median : 6.830   Median : 9.210   Median : 8.080   Median : 8.170  
 Mean   : 7.092   Mean   : 9.797   Mean   : 8.494   Mean   : 8.496  
 3rd Qu.: 9.670   3rd Qu.:12.960   3rd Qu.:11.420   3rd Qu.:11.210  
 Max.   :26.160   Max.   :30.370   Max.   :31.080   Max.   :25.880  
      CLO              BEL             MAL       
 Min.   : 0.040   Min.   : 0.13   Min.   : 0.67  
 1st Qu.: 5.330   1st Qu.: 8.71   1st Qu.:10.71  
 Median : 8.290   Median :12.50   Median :15.00  
 Mean   : 8.707   Mean   :13.12   Mean   :15.60  
 3rd Qu.:11.630   3rd Qu.:16.88   3rd Qu.:19.83  
 Max.   :28.210   Max.   :42.38   Max.   :42.54  

wind> wind.loc
         Station Code       Latitude     Longitude MeanWind
1       Valentia  VAL        51d56'N       10d15'W     5.48
2      Belmullet  BEL        54d14'N       10d00'W     6.75
3    Claremorris  CLA        53d43'N        8d59'W     4.32
4        Shannon  SHA        52d42'N        8d55'W     5.38
5  Roche's Point  RPT        51d48'N        8d15'W     6.36
6           Birr  BIR        53d05'N        7d53'W     3.65
7      Mullingar  MUL        53d32'N        7d22'W     4.38
8     Malin Head  MAL        55d22'N        7d20'W     8.03
9       Kilkenny  KIL        52d40'N        7d16'W     3.25
10        Clones  CLO        54d11'N        7d14'W     4.48
11        Dublin  DUB        53d26'N        6d15'W     5.05
12       Roslare  ROS 52d16'56.791"N 6d21'25.056"W     6.00

wind> wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]])))

wind> wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]])))

wind> coordinates(wind.loc) = ~x+y

wind> # fig 1:
wind> ## No test: 
wind>  if (require(mapdata)) {
wind+    map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5))
wind+    plot(wind.loc, add=TRUE, pch=16)
wind+    text(coordinates(wind.loc), pos=1, label=wind.loc$Station)
wind+  }
Loading required package: mapdata

wind> ## End(No test)
wind> wind$time = ISOdate(wind$year+1900, wind$month, wind$day)

wind> # time series of e.g. Dublin data:
wind> plot(DUB~time, wind, type= 'l', ylab = "windspeed (knots)", main = "Dublin")
Hit <Return> to see next plot: 

wind> # fig 2:
wind> #wind = wind[!(wind$month == 2 & wind$day == 29),]
wind> wind$jday = as.numeric(format(wind$time, '%j'))

wind> windsqrt = sqrt(0.5148 * as.matrix(wind[4:15]))

wind> Jday = 1:366

wind> windsqrt = windsqrt - mean(windsqrt)

wind> daymeans = sapply(split(windsqrt, wind$jday), mean)

wind> plot(daymeans ~ Jday)
Hit <Return> to see next plot: 

wind> lines(lowess(daymeans ~ Jday, f = 0.1))

wind> # subtract the trend:
wind> meanwind = lowess(daymeans ~ Jday, f = 0.1)$y[wind$jday]

wind> velocity = apply(windsqrt, 2, function(x) { x - meanwind })

wind> # match order of columns in wind to Code in wind.loc:
wind> pts = coordinates(wind.loc[match(names(wind[4:15]), wind.loc$Code),])

wind> # fig 3, but not really yet...
wind> dists = spDists(pts, longlat=TRUE)

wind> corv = cor(velocity)

wind> sel = !(as.vector(dists) == 0)

wind> plot(as.vector(corv[sel]) ~ as.vector(dists[sel]),
wind+ 	xlim = c(0,500), ylim = c(.4, 1), xlab = "distance (km.)", 
wind+ 	ylab = "correlation") 
Hit <Return> to see next plot: 

wind> # plots all points twice, ignores zero distance 
wind> 
wind> # now really get fig 3:
wind> ros = rownames(corv) == "ROS"

wind> dists.nr = dists[!ros,!ros]

wind> corv.nr = corv[!ros,!ros]

wind> sel = !(as.vector(dists.nr) == 0)

wind> plot(as.vector(corv.nr[sel]) ~ as.vector(dists.nr[sel]), pch = 3,
wind+ 	xlim = c(0,500), ylim = c(.4, 1), xlab = "distance (km.)", 
wind+ 	ylab = "correlation") 
Hit <Return> to see next plot: 

wind> # add outlier:
wind> points(corv[ros,!ros] ~ dists[ros,!ros], pch=16, cex=.5)

wind> xdiscr = 1:500

wind> # add correlation model:
wind> lines(xdiscr, .968 * exp(- .00134 * xdiscr))

> m = map2SpatialLines(
+ 	map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5), plot=F))
Loading required package: maps
Error in cbind(map$x, map$y) : could not find function "map"
In addition: Warning messages:
1: In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘mapdata’
2: In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘maps’