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
116
117
118
119
120
121
122
| heure <- c("05:17:53","17:53:48","18:10:49","20:04:30","21:03:18","22:21:14")
df <- data_frame(
heure = as_hms(heure),
latitude =c(46.50253,46.24055,46.22687,46.22042,46.20115,46.16915),
longitude = c(-1.788917,-1.360150,-1.337583,-1.317250,-1.279983,-1.245900),
vitesse =c(5.54,7.23,2.28,3.92,7.23,5.20),
angle = c(231,290,332,223,283,220)
)
#Step 1 - Prepare data frame (`df`) with georeferenced data: wind speed and wind direction.
# Setting existing coordinate as lat-long system (WSG84)
cord.dec = SpatialPoints(cbind(df$longitude, df$latitude), proj4string=CRS("+proj=longlat"))
# Transforming coordinate to WSG84 UTM zone 29 using EPSG=2154 (france)
cord.UTM <- spTransform(cord.dec, CRS("+proj=utm +zone=29N ellps=WGS84"))
#Dataframe with georeferenced wind data
id <- c(1:length(df$angle))
df2 <- data.frame(id=id,
start.x= cord.UTM$coords.x1,
start.y= cord.UTM$coords.x2,
w.speed= df$vitesse,
w.direction= df$angle,
w.temps= df$heure)
#Step 2 - Complement `df` with auxiliary coordinates for representing wind as arrowhead lines.
#Line parameters
line.length <- 1000 #length of polylines representing wind in the map (meters)
arrow.length <- 300 #lenght of arrowhead leg (meters)
arrow.angle <- 120 #angle of arrowhead leg (degrees azimuth)
#Generate data frame with auxiliary coordinates
end.xy.df2 <- data.frame(end.x=NA,end.y=NA,end.arrow.x=NA,end.arrow.y=NA)
for (i in c(1:nrow(df2))){
#coordinates of end points for wind lines (the initial points are the ones where data was observed)
if (df2$w.direction[i] <= 90) {
end.x <- df2$start.x[i] + (cos((90 - df2$w.direction[i]) * 0.0174532925) * line.length)
} else if (df2$w.direction[i] > 90 & df2$w.direction[i] <= 180) {
end.x <- df2$start.x[i] + (cos((df2$w.direction[i] - 90) * 0.0174532925) * line.length)
} else if (df2$w.direction[i] > 180 & df2$w.direction[i] <= 270) {
end.x <- df2$start.x[i] - (cos((270 - df2$w.direction[i]) * 0.0174532925) * line.length)
} else {end.x <- df2$start.x[i] - (cos((df2$w.direction[i] - 270) * 0.0174532925) * line.length)}
if (df2$w.direction[i] <= 90) {
end.y <- df2$start.y[i] + (sin((90 - df2$w.direction[i]) * 0.0174532925) * line.length)
} else if (df2$w.direction[i] > 90 & df2$w.direction[i] <= 180) {
end.y <- df2$start.y[i] - (sin((df2$w.direction[i] - 90) * 0.0174532925) * line.length)
} else if (df2$w.direction[i] > 180 & df2$w.direction[i] <= 270) {
end.y <- df2$start.y[i] - (sin((270 - df2$w.direction[i]) * 0.0174532925) * line.length)
} else {end.y <- df2$start.y[i] + (sin((df2$w.direction[i] - 270) * 0.0174532925) * line.length)}
#coordinates of end points for arrowhead leg lines (the initial points are the previous end points)
end.arrow.x <- end.x + (cos((df2$w.direction[i] + arrow.angle) * 0.0174532925) * arrow.length)
end.arrow.y <- end.y - (sin((df2$w.direction[i] + arrow.angle) * 0.0174532925) * arrow.length)
end.xy.df2 <- rbind(end.xy.df2,c(end.x,end.y,end.arrow.x,end.arrow.y))
}
end.xy2 <- end.xy.df2[-1,]
df2 <- data.frame(df2,end.xy2) #df with observed and auxiliary variables
#Step 3 - Create an object of class `SpatialLinesDataFrame` to use within `leaflet`.
lines <- data.frame(cbind(lng=c(df2$start.x,df2$end.x,df2$end.arrow.x),
lat=c(df2$start.y,df2$end.y,df2$end.arrow.y),
id=c(rep(df2$id,3))))
lines.list <- list()
for (i in c(1:max(lines$id))){
line <- subset(lines,lines$id==i)
line <- as.matrix(line[,c(1:2)])
line <- Line(line) #object of class 'Line'
lines.list[[i]] <- Lines(list(line), ID = i) #list of 'objects'Lines'
}
sp.lines <- SpatialLines(lines.list) #object of class 'SpatialLines'
proj4string(sp.lines) <- CRS("+init=epsg:3857")
#Convert CRS to geographic coordinates (http://spatialreference.org/ref/epsg/4326/)
#for overlaying on OpenStreetMaps tiles in Leaflet
sp.lines <- spTransform(sp.lines, CRS("+init=epsg:4326"))
rownames(df2) = df2$id
#Join wind variables (id, speed, direction and date) to object of class 'SpatialLines'
sp.lines.df <- SpatialLinesDataFrame(sp.lines, df2[,c(1,4:6)]) #object of class 'SpatialLinesDataFrame'
str(sp.lines.df) #inspect object structure
#Step 4 - Generate interactive and **static** map of wind speed and direction.
#popup settings
labels <- paste0("Vitesse vent réel: ",sp.lines.df@data$w.speed," knots/h<br>",
"Angle vent réel : ",sp.lines.df@data$w.direction," degrees<br>",
"heure: ", sp.lines.df@data$w.temps)
#pallete settings
pal <- colorNumeric(palette = colorRampPalette(c("red", "blue"))(5),
domain = 0:10)
#Create object fo class 'leaflet' 'htmlwidget'
m <- leaflet(sp.lines.df) %>%
addTiles() %>% # add default OpenStreetMap map tiles
addPolylines(color = ~pal(w.speed), opacity=1, weigh = 3, popup = labels) %>%
addLegend("bottomright", pal = pal, values = ~w.speed,
title = "Vitesse vent réel <br> (knots/h)",
opacity = 1) %>%
fitBounds(sp.lines.df@bbox[1,1], sp.lines.df@bbox[2,1], sp.lines.df@bbox[1,2], sp.lines.df@bbox[2,2])
#Plot map
m |
mauvais référentiel de projection géographique sous R
Répondre avec citation
Partager