first_glance

Scientific Questions

Routing algorithms in the ocean are hard (Christiansen, Fagerholt, & Ronen, 2004).

What kind of questions and results are we looking for?

• Q: How can we improve routing from a simple grid?

• A: [Generic] Use triangulated network. Compare distances between start/end points. Also compare computational time.

• A: [Application] Compare with extra whale risk conserved & less shipping distance traversed for this example. Results:

  • Show boxplot of grid vs triangulated routing and efficiencies gained on y axis for whales & shipping

  • Show plot of how parameter for triangulation (like angle) on x-axis changes, how efficincies are gained in whale conservation / shipping distance on y-axis

  • Showing the Shiny app before/after this new routing mechanism. Embedded as a HTML slideshow.

• Q: How can we get faster results for the shortest path problem than Dijkstra’s algorithm?

• A: [Generic] Use the A* algorithm, or one of its variants depending on how large the graph is. We can create large sample graphs to compare times and solutions for different-sized graphs.

• A: [Generic] Look into UAV path-planning methodologies, since many of those applications have continuous paths rather than discrete.

Technical Questions

Let’s try to use https://github.com/dgrimsman/boot-camp/issues.

Data

British Columbia (BC)

• land: bc_spp_gcs.shp
• whales: v72zw_epsg3857.grd

raster of whale hot spots

library(raster) # install.packages('rgdal')

whales = raster('data/v72zw_epsg3857.grd')
plot(whales)

ports to route between

library(readr)

ports = read_csv('data/ports_bc.csv')

interactive map

library(leaflet)

leaflet() %>%
  addProviderTiles("Stamen.TonerLite", options = providerTileOptions(noWrap = TRUE)) %>% 
  addRasterImage(
    whales, opacity = 0.8, project = F, group='Raster',
    colors = colorNumeric(
      palette = 'Reds', domain = c(cellStats(whales, 'min'), cellStats(whales, 'max')), 
      na.color = "#00000000", alpha = TRUE)) %>%
  addCircleMarkers(
    ~lon, ~lat, color='blue', data=ports, layerId=~name, group='Ports',
    popup = ~sprintf('<b>%s</b><br>%0.2f, %0.2f', name, lon, lat)) %>%
  addLayersControl(
    overlayGroups = c('Raster', 'Ports'),
    options = layersControlOptions(collapsed=T))

Testing New Packages

RBGL

source("https://bioconductor.org/biocLite.R")
biocLite("RBGL")

library(RBGL)

# astar search
con <- file(system.file("XML/dijkex.gxl",package="RBGL"), open="r")
coex <- fromGXL(con)
close(con)
astarSearch(coex)

RTriangle

• Computational Geometry in Python

• Triangle: .poly files

library(RTriangle)
library(rgdal)

## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
##  Path to GDAL shared files: /usr/share/gdal/1.10
##  Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-3

library(rgeos)

## rgeos version: 0.3-20, (SVN revision 535)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-3 
##  Polygon checking: TRUE

## 
## Attaching package: 'rgeos'

## The following object is masked from 'package:RTriangle':
## 
##     triangulate

library(RTriangle)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(Matrix)
library(igraph)

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:dplyr':
## 
##     %>%, as_data_frame, groups, union

## The following object is masked from 'package:rgeos':
## 
##     union

## The following object is masked from 'package:leaflet':
## 
##     %>%

## The following object is masked from 'package:raster':
## 
##     union

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

library(sp)
library(gdistance)

## 
## Attaching package: 'gdistance'

## The following object is masked from 'package:igraph':
## 
##     normalize

triangulate = RTriangle::triangulate # override rgeos::triangulate

## Create an object with a concavity
p <- pslg(P=rbind(c(0, 0), c(0, 1), c(0.5, 0.5), c(1, 1), c(1, 0)),
          S=rbind(c(1, 2), c(2, 3), c(3, 4), c(4, 5), c(5, 1)))

## Plot it
plot(p)

## Triangulate it
tp <- triangulate(p)
plot(tp)

## Load a data set containing a hole
A <- read.pslg(file.path(system.file(package = "RTriangle"), "extdata", "A.poly"))
plot(A)

## Triangulate the PSLG
tA <- triangulate(A)
plot(tA)

## Triangulate the PSLG with triangles in which no angle
## is smaller than 20 degrees
tA <- triangulate(A, q=20)
plot(tA)

## Triangulate the PSLG with triangles in which no triangle has 
## area greater than 0.001
tA <- triangulate(A, a=0.001)
plot(tA)

source('custom.R')
ply = read.pslg('data/superior.poly')
plot(ply)

t_ply <- triangulate(ply)
plot(t_ply)

t_ply <- triangulate(ply, q=20, D=T) # 'pq20D'
plot(t_ply)

Triangulate BC

Triangularize BC raster · Issue #1 · dgrimsman/boot-camp

Get Extent Polygon of Raster

whales_1 = whales
whales_1[!is.na(whales_1)] = 1
plot(whales_1)
whales_p = rasterToPolygons(whales_1, n=8, na.rm=T, digits=12, dissolve=T)
plot(whales_p)

writeOGR(whales)
writeOGR(whales_p, 'data/bc', layer='layer', driver='GeoJSON')
file.rename('data/bc','data/bc.geojson')

whales_p = readOGR(dsn='data/bc.geojson', layer='OGRGeoJSON')

#names(whales_p)
#slotNames(whales_p)
#whales_p@polygons

# vertices = data_frame(i=numeric(), x=numeric(), y=numeric())
# for (i in 1:length(whales_p@polygons[[1]]@Polygons)){ # p = whales_p@polygons[[1]]@Polygons[1]
#   vertices = bind_rows(
#     vertices,
#     a = whales_p@polygons[[1]]@Polygons[[i]]@coords %>%
#       as.data.frame() %>%
#       mutate(i = i))
# }

vertices = data_frame(i=numeric(), x=numeric(), y=numeric())
segments = data_frame(i=numeric(), start=numeric(), fin=numeric())
start_v = 1
fin_v = 1
for (i in 1:length(whales_p@polygons[[1]]@Polygons)){ # p = whales_p@polygons[[1]]@Polygons[1]
  p = whales_p@polygons[[1]]@Polygons[[i]]
  if (p@area < 1e10){
    next
  }
  
  # save vertices
  coords = rename(as.data.frame(p@coords), c('V1'='x', 'V2'='y'))
  vertices = bind_rows(vertices, coords)
  vertices = head(vertices, -1) # because last vertex is redundant
  
  #save segments
  num_vert = dim(p@coords)[1] - 1
  fin_v = start_v + num_vert - 1
  start = c(start_v:fin_v)
  fin = c((start_v + 1):fin_v, start_v)
  segments = bind_rows(
    segments,
    a = data_frame(start, fin))
  start_v = fin_v + 1
}
vertices$i<-seq.int(nrow(vertices))
segments$i<-seq.int(nrow(segments))

idx = which(duplicated(vertices[,c('x','y')]))
vertices$x[idx] = vertices$y[idx] + runif(length(idx), min=0.01, max=0.03)

# http://www.cs.cmu.edu/~quake/triangle.poly.html
# DONE: @bbest
# First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>
# Following lines: <vertex #> <x> <y> [attributes] [boundary marker]
poly = 'data/bc.poly'
write_file(sprintf('%d 2 0 0\n', nrow(vertices)), poly)
for (i in 1:nrow(vertices)){ # i =1
  row = vertices[i,]
  write_file(sprintf('%d %0.2f %0.2f\n', i, row$x, row$y), poly, append=T)
}

# DONE: @dgrimsman
# One line: <# of segments> <# of boundary markers (0 or 1)>
# Following lines: <segment #> <endpoint> <endpoint> [boundary marker]
write_file(sprintf('%d 0\n', nrow(segments)), poly, append=T)
for (i in 1:nrow(segments)){ # i =1
  row = segments[i,]
  write_file(sprintf('%d %g %g\n', i, row$start, row$fin), poly, append=T)
}

# TODO: @dgrimsman
# One line: <# of holes>
# Following lines: <hole #> <x> <y>
# # holes = data_frame(i=numeric(), x=numeric(), y=numeric())
# write_file(sprintf('%d\n', length(whales_p@polygons[[1]]@Polygons) - 1), poly, append=T)
# j = 1
# for (i in 1:length(whales_p@polygons[[1]]@Polygons)){
#   p = whales_p@polygons[[1]]@Polygons[[i]]
#   if (p@area <1e10){
#     hole_coord = p@labpt
#     x = hole_coord[1]
#     y = hole_coord[2]
#     # holes = rbind(holes, c(i, hole_coord[1, 1], hole_coord[1, 2]))
#     write_file(sprintf('%d %g %g\n', j, x, y), poly, append=T)
#     j = j + 1
#   }
# }

write_file(sprintf('%d\n', 0), poly, append=T)

bc_poly = read.pslg('data/bc.poly') # not yet working b/c segments needed
plot(bc_poly)

bc_poly_tri <- triangulate(bc_poly, q=20, D=T) # 'pq20D'
plot(bc_poly_tri)

# create adjacency matrix
num_v_tri = dim(bc_poly_tri$P)[1]
adj_mat <- Matrix(0, nrow = num_v_tri, ncol = num_v_tri, sparse=TRUE)
for (i in 1:dim(bc_poly_tri$E)[1]){
  edge = bc_poly_tri$E[i,]
  p1 = edge[1]
  p2 = edge[2]
  d = sqrt((bc_poly_tri$P[p1][1] - bc_poly_tri$P[p2][1])^2 + (bc_poly_tri$P[p1][2] - bc_poly_tri$P[p2][2])^2)
  adj_mat[p1, p2] = d
  adj_mat[p2, p1] = d
}

#create graph and find shortest path
graph = graph_from_adjacency_matrix(adj_mat, weighted=TRUE, mode="undirected")
path_tri = shortest_paths(graph, 1, 1000)

# ones map for getting linear path
r1 = whales
r1[!is.na(whales)] = 1

lonlat1 = SpatialPoints(data.frame(-127.08, 50.59))
lonlat2 = SpatialPoints(data.frame(-130.32, 54.31))

xy = c(coordinates(lonlat1), coordinates(lonlat2)) %>%
      matrix(ncol=2, byrow=T) %>%
      SpatialPoints(crs("+init=epsg:4326")) %>%
      spTransform(crs("+init=epsg:3857"))

# update to nearest non-NA points on raster for shortestPath to work
xy = c(
  coordinates(r1)[which.min(mask(distanceFromPoints(r1, xy[1]), r1)),],
  coordinates(r1)[which.min(mask(distanceFromPoints(r1, xy[2]), r1)),]) %>%
  matrix(ncol=2, byrow=T) %>%
  SpatialPoints(crs("+init=epsg:3857"))

path_grid = shortestPath(
    #geoCorrection(transition(1 / r1, mean, directions=8)), 
    geoCorrection(transition(1 / (r1), mean, directions=8), type="c"), 
    xy[1],
    xy[2],
    output='SpatialLines')
plot(path_grid)

R Tips Background

• HTML Documents - rmarkdown
• Bibliographies and Citations
• bbest/rmarkdown-example: R Markdown example showing figures & tables with captions, equations, inline R values and references with a Zotero library

References

Christiansen, M., Fagerholt, K., & Ronen, D. (2004). Ship routing and scheduling: Status and perspectives. Transportation Science, 38(1), 1–18. Retrieved from http://pubsonline.informs.org/doi/abs/10.1287/trsc.1030.0036