-
Notifications
You must be signed in to change notification settings - Fork 0
/
GlobCover_squirrel.Rmd
187 lines (143 loc) · 6.78 KB
/
GlobCover_squirrel.Rmd
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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
---
title: "Shortest path across Spain: the squirrel problem"
output:
html_notebook:
toc: yes
toc_collapse: no
toc_float: yes
html_document:
df_print: paged
toc: yes
---
```{r, message=FALSE, warning=FALSE, include=FALSE}
# Load packages
library(raster)
library(leaflet)
library(plotly)
library(igraph)
library(Matrix)
library(htmlwidgets)
# Load homemade functions
source('build_matrix_graph.R')
```
This notebook gathers the necessary analyses for the [squirrel dashboard](http://jsaezgallego.com/GlobCover_maps_squirrel/). The notebook consists of the following steps:
1. Download data from the internet. The data consists of GIS information on whether a piece of land is a forest, a river, a man-made construction, etc.
2. Create path matrix. The squirrel takes one step at a time. In terms of a raster map, this means we can reach each pixel only from adjacent pixels. Think of it as a huge spare matrix.
3. Optimize path. Calculate shortest path from A to B. The cost of each segment is the pixel value (forest = no cost to move)
4. Visualize solution
# Download and load raster
The raster is freely available from the web
```{r, echo=TRUE}
# Download globCover dataset
# download.file('http://due.esrin.esa.int/files/Globcover2009_V2.3_Global_.zip',destfile = "Globcover2009_V2.3_Global_.zip")
# unzip('Globcover2009_V2.3_Global_.zip',exdir = 'GlobCover')
# Load raster
glob = raster("GlobCover/GLOBCOVER_L4_200901_200912_V2.3.tif")
```
We also need to map the Globcover labels to a continous number, as indicated by the legend found [here](http://due.esrin.esa.int/page_globcover.php). We set water bodies to _NA_. Later on, we set the roughness ot water bodies to 100 to indicate that the squirrel cannot swim over them.
```{r}
# map labels to roughness values
globCover = c(11,14,20,30,40,50,60,70,90,100,
110,120,130,140,150,160,170,180,190,200,210,220,230)
roughness = c(0.10, 0.10, 0.30, 0.30, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,
1.50, 0.50, 0.10, 0.03, 0.05, 0.50, 0.60, 0.20, 1.00, NA, NA, NA, NA)
```
## Map of Spain
We chop a small section of the world map around the Iberian peninsula. The code is made to work only for squared rasters.
The start and end points of the path are set to be the northernmost point (Estaca do bares) and the sourthernmost one (Tarifa).
```{r, message=FALSE, warning=FALSE}
# Chop a small section
lat = 40
lon = -5.1
n=4.5 # width of the box
cropbox2 <-c(lon-n,lon+n,lat-n,lat+n) # crop a square box
spain_img = crop(glob,cropbox2)
# re-classify the raster values to numeric
small_img_plot = reclassify(spain_img,cbind(globCover,roughness))
# Let's tranform NA values to "-100" roughness so that we squirrel cannot swim. This one is only used for calculations.
small_img = reclassify(small_img_plot,c(NA,-100,-100))
# Save the rasters for later use in the dashboard
raster::writeRaster(spain_img, "Files/GlobCover_Spain.tif",overwrite=TRUE)
# Plot
# plot(spain_img) # Original globcover labels
plot(small_img_plot) # transformed labels, with NA
# plot(small_img) # we use this one for the analysis with NA=-100
# Start and end points
start_point = c(43.780824, -7.681479) # Estaca do bares
end_point = c(36.003797, -5.609124)# Tarifa
# Plot start and end points
points(x=start_point[2],y=start_point[1],pch=19)
points(x=end_point[2],y=end_point[1],pch=19)
```
Dimensions of the raster:
```{r, echo=TRUE, message=FALSE, warning=FALSE}
# translate raster to matrix
mat = raster::as.matrix(small_img)
# create 0-1 edge matrix
N = dim(mat)[1] * dim(mat)[2] # Number of vertices
Ni= dim(mat)[1] # number of rows
Nj= dim(mat)[2] # number of columns
N
```
# Path matrix
Next we create the graph as a matrix. The matrix is equal to a 0 if two pixels are not linked, and equal to the roughness of theincoming pixel if they are. Two pixels are linked if they are adjacent to each other. This huge matrix of dimension N by N is really sparse and unless optimized takes too long to be created. It is build inside the function _build_matrix_graph()_ and it took few tries to make it run at reasonable speeds.
The values of the raster are slighly modified to be used in a standard shortest-path algorithm: 0.0000001 weight on an edge means "forest", while 1.5 means "bare land" and 100 means "water bodies".
```{r, echo=TRUE, message=FALSE, warning=FALSE}
# Call our homemade function
val_raster = 1.5-values(small_img)+0.0000001
t = Sys.time()
mat_path = build_matrix_graph(Ni,Nj,val_raster = val_raster)
Sys.time() - t
mat_path[1:10,1:10]
```
# Shortest path
The shortest path is claulcated using the [igraph](http://igraph.org/) library. It works surprisingly fast.
```{r, echo=TRUE, message=FALSE, warning=FALSE}
# Build igraph
g <- graph.adjacency(mat_path, mode="directed",weighted = T)
## Calculate shortest path between two poitns
# Find index of the start and end point
img_points = rasterToPoints(small_img)
idx_start = which.min( abs( img_points[,1] - start_point[2]) + abs( img_points[,2] - start_point[1]))
idx_end = which.min( abs( img_points[,1] - end_point[2]) + abs( img_points[,2] - end_point[1]))
# shortest path algorithm
news.path <- shortest_paths(g,from = idx_start,
to = idx_end,
output = "both") # both path nodes and edges
# Find those poitns in lat,lon
path_squirrel = news.path$vpath[[1]] %>% as.numeric
path_points = img_points[path_squirrel,1:2]
# Save squirrel path
write.csv(data.frame(lon = path_points[,1],
lat = path_points[,2],
n = path_squirrel),
file="Files/path_coordinates_solution.csv",row.names=F)
```
# Visualize solution
The path and the roughness (land cover) map is displayed below. Due to its size, the land cover map needs to be aggregated by a factor of 3(i.e., reduce its size).
```{r, echo=FALSE, message=FALSE, warning=FALSE}
# The raster is too big to be displayed so let's shrink it
agg_small_img_plot = raster::aggregate(small_img_plot,fact=3)
# Make a leaflet visualization
pal2 = colorNumeric("RdYlGn", domain = NULL)
ll_map = leaflet() %>%
# Base maps
addProviderTiles('Esri.WorldGrayCanvas',group="Esri Grey") %>%
addProviderTiles('Esri.WorldTopoMap',group="Esri Topo") %>%
addProviderTiles("Esri.WorldImagery", group = "Esri Image") %>%
# # Roughness info
addRasterImage(agg_small_img_plot, opacity = 0.5,colors=pal2, group="Roughness map") %>%
addLegend(position="bottomleft", pal = pal2, values = values(agg_small_img_plot),
title = "Roughness",
opacity = 1) %>%
# path info
addPolylines(lng=path_points[,1],lat=path_points[,2],group="Path") %>%
# Control groups
addLayersControl(
baseGroups = c("Esri Grey","Esri Image","Esri Topo"),
overlayGroups = c("Roughness map", "Path"),
options = layersControlOptions(collapsed = F)
)
ll_map
# saveWidget(ll_map,'example.html')
```