Maps

Set options and load packages

install.packages("rgdal")
install.packages("maptools")
install.packages("rworldmap")
install.packages("sp")
install.packages("maps")
install.packages("rworldmap")
install.packages("rgdal")
install.packages("rgeos")
install.packages("raster")
install.packages("plyr")
require(devtools)
install_github("riatelab/cartography")
install.packages("maptools")
install.packages("PBSmapping")
install.packages("ggmap")
install.packages("geosphere")
install.packages("RgoogleMaps")
install.packages("OpenStreetMap")
install.packages("tripack")

options(digits=3, width=70)
#install.packages("rgdal")
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.2-15, (SVN revision 691)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /usr/share/gdal/2.1
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5

library(maptools)

## Checking rgeos availability: TRUE

Proj <- CRS("+proj=longlat +datum=WGS84")
download.file("http://freakonometrics.free.fr/CAN1.zip",destfile="CAN.zip")
unzip("CAN.zip")
CAN_shp <- readShapePoly("CAN_adm1.shp", verbose=TRUE, proj4string=Proj)

## Warning: use rgdal::readOGR or sf::st_read

## Shapefile type: Polygon, (5), # of Shapes: 13

plot(CAN_shp)

new_CAN_shp <- spTransform(CAN_shp, CRS("+init=epsg:26978"))
plot(new_CAN_shp)

library(sp)
library(rworldmap)

## ### Welcome to rworldmap ###

## For a short introduction type :   vignette('rworldmap')

data(countriesLow)
crs.WGS84 <- CRS("+proj=longlat +datum=WGS84")
countries.WGS84 <- spTransform(countriesLow, crs.WGS84)
plot(countries.WGS84,col="light green")

crs.laea <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") 
countries.laea <- spTransform(countriesLow, crs.laea) 
plot(countries.laea,col="light green")

library(maps)
data("world.cities")
world.cities2 = world.cities[world.cities$pop>100000,]
plot(world.cities2$lon,world.cities2$lat,pch=19,cex=.7,axes=FALSE,xlab="",ylab="")

loc="https://raw.githubusercontent.com/jjrom/hydre/master/mapserver/geodata/rivers/GRDC_687_rivers"
ext=c(".dbf",".prj",".sbn",".sbx",".shp",".shp.xml",".shx")
for(i in 1:length(ext)){
 loc1=paste(loc,ext[i],"?raw=true",sep="")
 loc2=paste("GRDC_687_rivers",ext[i],sep="")
 download.file(loc1,loc2)
}
library(sp)
shape <- readShapeLines("GRDC_687_rivers.shp")

## Warning: use rgdal::readOGR or sf::st_read

plot(shape,col="blue")

#crs.laea <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") 
#shape.nouveau = spTransform(shape,crs.laea)
#plot(shape.nouveau,col="blue")

#download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_roads.zip", "routes.zip")
unzip("routes.zip")

## Warning in unzip("routes.zip"): erreur 1 lors de l'extraction d'un
## fichier zip

shape <- readShapeLines("ne_10m_roads.shp")

## Warning: use rgdal::readOGR or sf::st_read

plot(shape,pch=19,lwd=.7)

download.file("http://freakonometrics.free.fr/Italy.zip",destfile="Italy.zip")
unzip("Italy.zip")
library(rgdal)
ita1<-readOGR(dsn="/home/arthur/Italy/",layer="ITA_adm1",verbose=FALSE)
plot(ita1,col="light green")

ita1@data[,"NAME_1"][1:6]

## [1] Abruzzo        Apulia         Basilicata     Calabria      
## [5] Campania       Emilia-Romagna
## 20 Levels: Abruzzo Apulia Basilicata Calabria ... Veneto

pos<-which(ita1@data[,"NAME_1"]=="Sicily")
Poly_Sicily<-ita1[pos,]
plot(ita1,col="light green")
plot(Poly_Sicily,col="yellow",add=TRUE)

plot(ita1,col="light green")
plot(Poly_Sicily,col="yellow",add=TRUE)
abline(v=5:20,col="light blue")
abline(h=35:50,col="light blue")
axis(1)
axis(2)

pos<-which(ita1@data[,"NAME_1"] %in% c("Abruzzo","Apulia","Basilicata","Calabria","Campania","Lazio","Molise","Sardegna","Sicily"))
ita_south <- ita1[pos,]
ita_north <- ita1[-pos,]
plot(ita1)
plot(ita_south,col="red",add=TRUE)
plot(ita_north,col="blue",add=TRUE)

library(rgeos)

## rgeos version: 0.3-26, (SVN revision 560)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE

ita_s<-gUnionCascaded(ita_south)
ita_n<-gUnionCascaded(ita_north)
plot(ita1)
plot(ita_s,col="red",add=TRUE)
plot(ita_n,col="blue",add=TRUE)

plot(ita1,col="light green")
plot(Poly_Sicily,col="yellow",add=TRUE)
gCentroid(Poly_Sicily,byid=TRUE)

## SpatialPoints:
##       x    y
## 14 14.1 37.6
## Coordinate Reference System (CRS) arguments: +proj=longlat
## +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

points(gCentroid(Poly_Sicily,byid=TRUE),pch=19,col="red")

G <- as.data.frame(gCentroid(ita1,byid=TRUE))
plot(ita1,col="light green")
text(G$x,G$y,1:20)

c1 <-G[c(17,20,6,16,8,12,2),]
c2 <- G[c(13,11,1,5,3,4),]
L1<-SpatialLines(list(Lines(list(Line(c1)),"Line1")))
L2<-SpatialLines(list(Lines(list(Line(c2)),"Line2")))
L1@proj4string<-ita1@proj4string
L2@proj4string<-ita1@proj4string
cross_road <-gIntersection(L1,L2)
cross_road@coords

##      x    y
## 1 11.1 44.0
## 1 14.2 41.7

plot(ita1,col="light green")
plot(L1,col="red",lwd=2,add=TRUE)
plot(L2,col="blue",lwd=2,add=TRUE)
plot(cross_road,pch=19,cex=1.5,add=TRUE)

library(sp)
library(raster)
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/GBR_adm2.rds","GBR_adm2.rds")
UK=readRDS("GBR_adm2.rds")
UK@data[159,"HASC_2"]="GB.NR"
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/IRL_adm0.rds","IRL_adm0.rds")
IRL=readRDS("IRL_adm0.rds")
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")
plot(IRL,add=TRUE)

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds")
FR=readRDS("FRA_adm0.rds")
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")
plot(IRL,add=TRUE)
plot(FR,add=TRUE)

loc="https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/12/EU-referendum-result-data.csv"
referendum=read.csv(loc,header=TRUE,dec=".",sep=",",stringsAsFactors = FALSE)
referendum=referendum[c(3,6,13,14)]
library(plyr)

## 
## Attaching package: 'plyr'

## The following object is masked from 'package:maps':
## 
##     ozone

referendum=ddply(referendum,.(Region,HASC_code),summarise,Remain=sum(Remain),Leave=sum(Leave))
sum(referendum$Leave)/(sum(referendum$Leave)+sum(referendum$Remain))

## [1] 0.519

referendum=referendum[c(referendum$Region!="Northern Ireland"),]
referendum=referendum[c(referendum$HASC_code!="Gibraltar"),]
row.names(referendum)=seq(1,nrow(referendum),1)
leave_or_remain=cbind(referendum,"Brexit?"=0)
leave_or_remain[,"Brexit?"]=ifelse(leave_or_remain$Remain<leave_or_remain$Leave,rgb(1,0,0,.7),rgb(0,0,1,.7))
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Brexit"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Brexit"]="blue"}else{
map_data[i,"Brexit"]=as.character(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],'Brexit?'])
}
}
plot(UK, col = map_data$Brexit, border = "gray1", xlim = c(-4,-2), ylim = c(50, 59), main="How the UK has voted?", bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)
legend(-1,59,c("Leave","Remain"),fill=c("red","blue"),bty="n")

#require(devtools)
#install_github("riatelab/cartography")
library(cartography)

## Loading required package: sf

## Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3

cols <- carto.pal(pal1 = "red.pal",n1 = 5,pal2="green.pal",n2=5)
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Percentage_Remain"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Percentage_Remain"]=55.78}else{         map_data[i,"Percentage_Remain"]=100*leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]/(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]+leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Leave"])
}
}
plot(UK, col = "grey", border = "gray1", xlim = c(-4,-2), ylim = c(50, 59),bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)
choroLayer(spdf = UK,
df = map_data,
var = "Percentage_Remain",
breaks = seq(0,100,10),
col = cols,
legend.pos = "topright",
legend.title.txt = "",
legend.values.rnd = 2,
add = TRUE)

library (maptools)
library (rgdal)
library (PBSmapping)

## 
## -----------------------------------------------------------
## PBS Mapping 2.70.4 -- Copyright (C) 2003-2017 Fisheries and Oceans Canada
## 
## PBS Mapping comes with ABSOLUTELY NO WARRANTY;
## for details see the file COPYING.
## This is free software, and you are welcome to redistribute
## it under certain conditions, as outlined in the above file.
## 
## A complete user guide 'PBSmapping-UG.pdf' is located at 
## /home/arthur/R/x86_64-pc-linux-gnu-library/3.3/PBSmapping/doc/PBSmapping-UG.pdf
## 
## Packaged on 2017-06-28
## Pacific Biological Station, Nanaimo
## 
## All available PBS packages can be found at
## https://github.com/pbs-software
## 
## To see demos, type '.PBSfigs()'.
## -----------------------------------------------------------

download.file("http://freakonometrics.free.fr/paris-cartelec.zip","paris-cartelec.zip")
unzip("paris-cartelec.zip")
paris <- readShapeSpatial ( "paris-cartelec.shp" )

## Warning: use rgdal::readOGR or sf::st_read

## Warning: use rgdal::readOGR or sf::st_read

## Warning in Polygon(coords = crds): less than 4 coordinates in polygon

## Warning in Polygon(coords = crds): less than 4 coordinates in polygon

## Warning in Polygon(coords = crds): less than 4 coordinates in polygon

## Warning in Polygon(coords = crds): less than 4 coordinates in polygon

## Warning in Polygon(coords = crds): less than 4 coordinates in polygon

proj4string ( paris ) <- CRS( "+init=epsg:2154" )
paris <- spTransform(paris , CRS( "+proj=longlat +ellps=GRS80" ) )
plot ( paris )

length(paris)

## [1] 850

poly_paris <- SpatialPolygons2PolySet(paris)
point=c(2.33,48.849)
plot(point,xlim=point[1]+c(-.01,.01),ylim=point[2]+c(-.01,.01))
plot(paris,add=TRUE)

plot(point,xlim=point[1]+c(-.01,.01),ylim=point[2]+c(-.01,.01))
plot(paris,add=TRUE)
library(rgeos)
library(geosphere)
idx=unique(poly_paris$PID)
for(i in idx){
ct=centroid(as.matrix(poly_paris[poly_paris$PID==i,c("X","Y")]))
text(ct[1],ct[2],i,col="blue")}

which.poly <- function ( point ) {
 idx=unique(poly_paris$PID)
 idx_valid=NULL
 for ( i in idx ) {
 pip=point.in.polygon (point[ 1 ] , point[ 2 ] ,
 poly_paris [ poly_paris$PID==i , "X" ] ,
 poly_paris [ poly_paris$PID==i , "Y" ] )
 if ( pip >0) idx_valid=c ( idx_valid , i ) }
 return ( idx_valid ) }
which.poly(point)

## [1] 86

library(ggmap)

## Loading required package: ggplot2

library(geosphere)
#(Hotel <- geocode("8 sinaran drive singapore"))
Hotel = c(103.845,1.320279)
#(Airport <- geocode("singapore airport"))
Airport = c(103.9915,1.36442)
distHaversine(Airport,Hotel, r=6378.137)

## [1] 17

#mapdist(as.numeric(Hotel),as.numeric(Airport))

Loc=data.frame(rbind(Airport, Hotel))
library(ggmap)
library(RgoogleMaps)
CenterOfMap <- (apply(Loc,2,mean))
names(CenterOfMap)=c("lon","lat")
SG <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "terrain", source = "google")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

SGMap <- ggmap(SG)
SGMap 

SGMap <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "roadmap")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false

ggmap(SGMap)

SGMap <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "satellite")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false

ggmap(SGMap)

SGMap <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "toner", source = "stamen")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=terrain&sensor=false

## Map from URL : http://tile.stamen.com/toner/11/1613/1015.png

## Map from URL : http://tile.stamen.com/toner/11/1614/1015.png

## Map from URL : http://tile.stamen.com/toner/11/1615/1015.png

## Map from URL : http://tile.stamen.com/toner/11/1616/1015.png

## Map from URL : http://tile.stamen.com/toner/11/1613/1016.png

## Map from URL : http://tile.stamen.com/toner/11/1614/1016.png

## Map from URL : http://tile.stamen.com/toner/11/1615/1016.png

## Map from URL : http://tile.stamen.com/toner/11/1616/1016.png

## Map from URL : http://tile.stamen.com/toner/11/1613/1017.png

## Map from URL : http://tile.stamen.com/toner/11/1614/1017.png

## Map from URL : http://tile.stamen.com/toner/11/1615/1017.png

## Map from URL : http://tile.stamen.com/toner/11/1616/1017.png

ggmap(SGMap)

library(OpenStreetMap)
map <- openmap(c(lat= 51.516,   lon= -.141), c(lat= 51.511,   lon= -.133))
map <- openproj(map, projection = "+init=epsg:27700") 
plot(map)

library(OpenStreetMap,verbose=FALSE,quietly=TRUE,warn.conflicts=FALSE)
map <- openmap(c(lat= 51.516, lon= -.141), c(lat= 51.511, lon= -.133))
map <- openproj(map, projection = "+init=epsg:27700") 
plot(map)

loc="http://freakonometrics.free.fr/cholera.zip"
download.file(loc,destfile="cholera.zip")
unzip("cholera.zip", exdir="./")
library(maptools,verbose=FALSE,quietly=TRUE,warn.conflicts=FALSE)
deaths<-readShapePoints("Cholera_Deaths")

## Warning: use rgdal::readOGR or sf::st_read

head(deaths@coords)

##   coords.x1 coords.x2
## 0    529309    181031
## 1    529312    181025
## 2    529314    181020
## 3    529317    181014
## 4    529321    181008
## 5    529337    181006

plot(map)
points(deaths@coords,col="red", pch=19, cex=.7 )

X = deaths@coords
library(KernSmooth)

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

kde2d = bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])))

## Warning in bw.ucv(X[, 1]): minimum occurred at one end of the range

library(grDevices,verbose=FALSE,quietly=TRUE,warn.conflicts=FALSE)
clrs = colorRampPalette(c(rgb(0,0,1,0), rgb(0,0,1,1)), alpha = TRUE)(20)
plot(map)
points(deaths@coords,col="red", pch=19, cex=.7 )
image(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE,col=clrs)
contour(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE)

library(ggmap,verbose=FALSE,quietly=TRUE,warn.conflicts=FALSE)
get_london = get_map(c(-.137,51.513), zoom=17)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=51.513,-0.137&zoom=17&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

london = ggmap(get_london)
london

df_deaths <- data.frame(X)
library(sp,verbose=FALSE,quietly=TRUE,warn.conflicts=FALSE)
library(rgdal,verbose=FALSE,quietly=TRUE,warn.conflicts=FALSE)
coordinates(df_deaths)=~coords.x1+coords.x2
proj4string(df_deaths)=CRS("+init=epsg:27700") 
df_deaths = spTransform(df_deaths,CRS("+proj=longlat +datum=WGS84"))
london + geom_point(aes(x=coords.x1,y=coords.x2),data=data.frame(df_deaths@coords),col="red")

## Warning: Removed 9 rows containing missing values (geom_point).

Consider a collection of roads

dd <- data.frame(Id=c(1,1,1,2,2,2,3,3,4,4),
                 X=c(3,5,8,5,7,8,3,5,5,6),
                 Y=c(9,8,3,4,7,4,8,9,7,8))
plot(dd$X,dd$Y)
for(i in unique(dd[,1])) lines(dd$X[dd$Id==i],dd$Y[dd$Id==i])

Let us create a shapefile from those segments

library(shapefiles)

## Loading required package: foreign

## 
## Attaching package: 'shapefiles'

## The following objects are masked from 'package:foreign':
## 
##     read.dbf, write.dbf

library(geosphere)
library(rgeos)
extend_shp <- function(dd){
  pb <- txtProgressBar(min = 0, max = 100, style = 3)
  newdd1=NULL
  newdd2=NULL
  cpte=0
  W=names(which(table(dd$Id)>1))
  dd=dd[dd$Id %in% as.numeric(W),]
  for(no in unique(dd$Id)){
   cpte=cpte+1
    if(floor(cpte %%(length(unique(dd$Id))/100)) == 0) setTxtProgressBar(pb, round(100*cpte/length(unique(dd$Id))))
    G=NULL
    DD=dd[dd$Id==no,2:3]
    p1=c(min(DD[,1]),min(DD[,2]),max(DD[,1]),max(DD[,2]))
    P1=cbind(c(p1[1],p1[1],p1[3],p1[3],p1[1]),c(p1[2],p1[4],p1[4],p1[2],p1[2]))
    box_i=function(i){
      b=dd[dd$Id==i,2:3]
      p2=c(min(b[,1]),min(b[,2]),max(b[,1]),max(b[,2]))
      P2=cbind(c(p2[1],p2[1],p2[3],p2[3]),c(p2[2],p2[4],p2[4],p2[2]))
      ct=function(j) point.in.polygon(P2[j,1],P2[j,2],P1[,1],P2[,2])
      sum(Vectorize(ct)(1:4))
    }
    vu = unique(dd$Id)[-which(unique(dd$Id)==no)]
    vu=vu[which(Vectorize(box_i)(vu)>1)]
    for(t in 2:nrow(DD)){
      P=DD[t-(1:0),]
      L1=SpatialLines(list(Lines(list(Line(P)),"Line1")))
      for(i in vu){
        L2=SpatialLines(list(Lines(list(Line(dd[dd$Id==i,2:3])),"Line2")))
        Inters=gIntersection(L1,L2)
        if(length(Inters)>0){
            if(class(Inters)=="SpatialPoints") G=rbind(G,c(Inters@coords,t))
            if(class(Inters)=="SpatialLines") G=rbind(G,c(Inters@lines[[1]]@Lines[[1]]@coords,t))
        
      }}
    P2=DD[1,]
    for(t in 2:nrow(DD)){
      if(sum(G[,3]==t)==0) P2=rbind(P2,DD[t,])
      if(sum(G[,3]==t)==1) P2=rbind(P2,G[G[,3]==t,1:2],DD[t,])
      if(sum(G[,3]==t)>1){
        sG=G[G[,3]==t,1:2]
        if(DD[t-1,1]<DD[t,1]) sG=sG[order(sG[,1]),] 
        if(DD[t-1,1]>DD[t,1]) sG=sG[-order(sG[,1]),] 
        colnames(sG)=c("X","Y")
        P2=rbind(P2,sG,DD[t,])
      }
    }
    dd2=data.frame(Id=no,X=P2$X,Y=P2$Y)
    dd2=dd2[!duplicated(dd2),]
    newdd1=rbind(newdd1,dd2)
    dd2=data.frame(Id=no,X1=P2$X[1:(nrow(dd2)-1)],Y1=P2$Y[1:(nrow(dd2)-1)],
                         X2=P2$X[2:(nrow(dd2))],Y2=P2$Y[2:(nrow(dd2))])
    dd2=dd2[!duplicated(dd2),]
    newdd2=rbind(newdd2,dd2)
  }}
  LIST_KNOTS = newdd1[!duplicated(newdd1[,2:3]),2:3]
  LIST_KNOTS=data.frame(LIST_KNOTS,NO=1:nrow(LIST_KNOTS))
  newdd2$Tr = 1:nrow(newdd2)
  depart  = data.frame(Tr=newdd2[,"Tr"],Id=newdd2[,"Id"],X=newdd2[,"X1"],Y=newdd2[,"Y1"])
  arrivee = data.frame(Tr=newdd2[,"Tr"],Id=newdd2[,"Id"],X=newdd2[,"X2"],Y=newdd2[,"Y2"])
  depart =  merge(depart,LIST_KNOTS,all.x=TRUE)
  arrivee = merge(arrivee,LIST_KNOTS,all.x=TRUE)
  names(depart)[names(depart)=="NO"]="knot1"
  names(arrivee)[names(arrivee)=="NO"]="knot2"
  base=merge(newdd2,depart[,c("Tr","Id","knot1")],all.x=TRUE,by="Tr")
  base=merge(base,arrivee[,c("Tr","Id","knot2")],all.x=TRUE,by="Tr") 
  base$distance = distHaversine( base[,c("X1","Y1")] ,  base[,c("X2","Y2")] , r =6378.137)
  base2=base[,c("knot1","knot2","distance","Id")]
  nb=base2[,c("knot2","knot1","distance","Id")];names(nb)=c("knot1","knot2","distance","Id")
  base2=rbind(base2,nb)
  return(list(shp=newdd1,
              knots=LIST_KNOTS,
              connected=base2))}
dd2 <- extend_shp(dd)

## 
  |                                                                  
  |                                                            |   0%
  |                                                                  
  |===============                                             |  25%
  |                                                                  
  |==============================                              |  50%
  |                                                                  
  |=============================================               |  75%
  |                                                                  
  |============================================================| 100%

dd2

## $shp
##    Id    X    Y
## 1   1 3.00 9.00
## 2   1 4.00 8.50
## 3   1 5.00 8.00
## 4   1 8.00 3.00
## 5   1 3.00 9.00
## 6   1 4.00 8.50
## 7   1 5.00 8.00
## 8   1 5.38 7.38
## 9   1 6.26 5.89
## 10  1 8.00 3.00
## 11  2 5.00 4.00
## 12  2 6.26 5.89
## 13  2 7.00 7.00
## 14  2 8.00 4.00
## 15  2 5.00 4.00
## 16  2 6.26 5.89
## 17  2 7.00 7.00
## 18  2 8.00 4.00
## 19  3 3.00 8.00
## 20  3 4.00 8.50
## 21  3 5.00 9.00
## 22  4 5.00 7.00
## 23  4 6.00 8.00
## 
## $knots
##       X    Y NO
## 1  3.00 9.00  1
## 2  4.00 8.50  2
## 3  5.00 8.00  3
## 4  8.00 3.00  4
## 8  5.38 7.38  5
## 9  6.26 5.89  6
## 11 5.00 4.00  7
## 13 7.00 7.00  8
## 14 8.00 4.00  9
## 19 3.00 8.00 10
## 21 5.00 9.00 11
## 22 5.00 7.00 12
## 23 6.00 8.00 13
## 
## $connected
##    knot1 knot2 distance Id
## 1      1     2    123.3  1
## 2      2     3    123.4  1
## 3      3     4    648.3  1
## 4      1     2    123.3  1
## 5      2     3    123.4  1
## 6      3     5     80.9  1
## 7      5     6    191.8  1
## 8      6     4    375.5  1
## 9      7     6    253.2  2
## 10     6     8    147.6  2
## 11     8     9    351.9  2
## 12     7     6    253.2  2
## 13     6     8    147.6  2
## 14     8     9    351.9  2
## 15    10     2    123.4  3
## 16     2    11    123.3  3
## 17    12    13    156.8  4
## 18     2     1    123.3  1
## 19     3     2    123.4  1
## 20     4     3    648.3  1
## 21     2     1    123.3  1
## 22     3     2    123.4  1
## 23     5     3     80.9  1
## 24     6     5    191.8  1
## 25     4     6    375.5  1
## 26     6     7    253.2  2
## 27     8     6    147.6  2
## 28     9     8    351.9  2
## 29     6     7    253.2  2
## 30     8     6    147.6  2
## 31     9     8    351.9  2
## 32     2    10    123.4  3
## 33    11     2    123.3  3
## 34    13    12    156.8  4

Intersections of roads are now added, as knots of the network

plot(dd$X,dd$Y,col="white")
for(i in unique(dd[,1])) lines(dd$X[dd$Id==i],dd$Y[dd$Id==i])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="blue")
points(dd$X,dd$Y,pch=19,col="red")

plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white")
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO)

plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="light blue")
clr <- c(rep("black",nrow(dd2$connected)-1),"red")
for(i in 1:nrow(dd2$connected)) arrows(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
                                       dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]],
                                       length=.1,col=clr[i])

library(RColorBrewer)
darkcols <- brewer.pal(8, "Dark2")
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="light blue",xlab="",ylab="",axes=FALSE)
for(i in 1:nrow(dd2$connected)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
                                         dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]],
                                         col=darkcols[dd2[["connected"]]$Id[i]])
points(dd$X,dd$Y,pch=19,col=darkcols[dd$Id],cex=3)
text(dd$X,dd$Y,1:nrow(dd),col="white")

We can also visualize all knots of the network

plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in 1:(nrow(dd2$connected)/2)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
                                             dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]],
                                             col=darkcols[dd2[["connected"]]$Id[i]])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="black",cex=3)
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO,col="white")

Adding Agents

agents=data.frame(X=c(5.5,3.5,6.5,6.5,4.5),
                  Y=c(4,9,8,5.5,9),
                  Type=as.factor(c("Demand","Demand","Supply","Supply","Supply")))
forme=c(15,19,0,1)
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in 1:(nrow(dd2$connected)/2)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
                                             dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="white",cex=3)
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=1,col="black",cex=3)
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO,col="black")
points(agents$X,agents$Y,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3)

Minimal Distance (Euclidean)

plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in which(agents$Type=="Demand")) segments(agents[i,"X"],agents[i,"Y"],agents[agents$Type=="Supply","X"],agents[agents$Type=="Supply","Y"])
points(agents$X,agents$Y,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3)

Closest distance Demand-Supply

close_supply <- function(coord,knots=agents[agents$Type=="Supply",]){
  knots[which.min(distHaversine( coord[,c("X","Y")] ,  knots[,c("X","Y")] , r =6378.137)),]
}
list_demand <- which(agents$Type=="Demand")
close_supply(agents[list_demand[1],])

##     X   Y   Type
## 4 6.5 5.5 Supply

plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in which(agents$Type=="Demand")) segments(agents[i,"X"],agents[i,"Y"],agents[agents$Type=="Supply","X"],agents[agents$Type=="Supply","Y"])
for(i in list_demand){
     supply=unlist(close_supply(agents[i,]))
     segments(agents[i,"X"],agents[i,"Y"],supply["X"],supply["Y"],col=darkcols[3],lwd=6)
}
points(agents$X,agents$Y,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3.2)
text(agents$X,agents$Y,1:nrow(agents),pch=forme[agents$Type],col="white")

plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE)
text(.125,.5,"DEMAND",srt=90,col=darkcols[1])
text(.875,.5,"SUPPLY",srt=270,col=darkcols[2])
i=which(agents$Type=="Demand")
y1=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
i=which(agents$Type=="Supply")
y2=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
id=which(agents$Type=="Demand")
for(i in 1:sum(agents$Type=="Demand")){
  d=distHaversine( agents[id[i],c("X","Y")] ,  agents[agents$Type=="Supply",c("X","Y")] , r =6378.137)
  segments(.25,y1[i],.75,y2,col=c("black",darkcols[3])[1+(d==min(d))],lwd=c(1,6)[1+(d==min(d))])
  text(.5,(y2+y1[i])/2,round(d),pos=1,col=c("black",darkcols[3])[1+(d==min(d))])
}
i=which(agents$Type=="Demand")
points(rep(.25,length(i)),y1,
       pch=forme[1],col=darkcols[1],cex=3.2)
text(rep(.25,length(i)),y1,i,col="white")
i=which(agents$Type=="Supply")
points(rep(.75,length(i)),y2,
       pch=forme[2],col=darkcols[2],cex=3.2)
text(rep(.75,length(i)),y2,i,col="white")

Agents Location and Knots

Moving agents to knots of the network

close_knot=function(coord,knots=dd2$knots){
  knots[which.min(distHaversine( coord[,c("X","Y")] ,  knots[,c("X","Y")] , r =6378.137)),]
}
shift_agents=agents
shift_agents$Xk=shift_agents$X
shift_agents$Yk=shift_agents$Y
shift_agents$NO=0
for(i in 1:nrow(shift_agents)) shift_agents[i,c("Xk","Yk","NO")]=close_knot(agents[i,])
shift_agents

##     X   Y   Type   Xk   Yk NO
## 1 5.5 4.0 Demand 5.00 4.00  7
## 2 3.5 9.0 Demand 3.00 9.00  1
## 3 6.5 8.0 Supply 6.00 8.00 13
## 4 6.5 5.5 Supply 6.26 5.89  6
## 5 4.5 9.0 Supply 5.00 9.00 11

plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in 1:(nrow(dd2$connected)/2)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
                                             dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="white",cex=3)
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=1,col="black",cex=3)
points(agents$Xk,agents$Yk,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3)
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO,col="black")
points(shift_agents$Xk,shift_agents$Yk,pch=forme[shift_agents$Type],col=darkcols[shift_agents$Type],cex=3.2)
points(shift_agents$X,shift_agents$Y,pch=forme[2+as.numeric(shift_agents$Type)],col=darkcols[as.numeric(shift_agents$Type)],cex=3)

Minimal Distance on a Network

plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE)
text(.125,.5,"DEMAND",srt=90,col=darkcols[1])
text(.875,.5,"SUPPLY",srt=270,col=darkcols[2])
i=which(agents$Type=="Demand")
y1=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
i=which(agents$Type=="Supply")
y2=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
id=which(agents$Type=="Demand")
for(i in 1:sum(agents$Type=="Demand")){
  d=distHaversine( agents[id[i],c("X","Y")] ,  agents[agents$Type=="Supply",c("X","Y")] , r =6378.137)
  segments(.25,y1[i],.75,y2,col=c("black",darkcols[3])[1+(d==min(d))],lwd=c(1,6)[1+(d==min(d))])
  text(.5,(y2+y1[i])/2,round(d),pos=1,col=c("black",darkcols[3])[1+(d==min(d))])
}
i=which(agents$Type=="Demand")
points(rep(.25,length(i)),y1,
       pch=forme[1],col=darkcols[1],cex=3.2)
text(rep(.25,length(i)),y1,i,col="white")
i=which(agents$Type=="Supply")
points(rep(.75,length(i)),y2,
       pch=forme[2],col=darkcols[2],cex=3.2)
text(rep(.75,length(i)),y2,i,col="white")

Application on a real network : subway in Paris

download.file("http://freakonometrics.free.fr/arcs.csv","arcs.csv")
arcs = read.csv(file = "arcs.csv",sep=";", header=FALSE) 
download.file("http://freakonometrics.free.fr/nodes.csv","nodes.csv")
nodes = read.csv("nodes.csv",sep=";", header=FALSE)
plot(nodes[1:302,c("V3","V4")])

subarcs=arcs[(arcs$V1<303)&(arcs$V2<303),]
plot(nodes[1:302,c("V3","V4")],pch=19,xlab="",ylab="",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
                                nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"])

library(OpenStreetMap)
map=openmap(c(lat=48.9,lon=2.26),c(lat=48.8,lon=2.45))
map=openproj(map,projection="+init=epsg:4326")
plot(map)
points(nodes[1:302,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
                                nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"])

Consider for instance locations of some hospitals in Paris,

library(ggmap)
adresses <- c("47 boulevard de l'hopital Paris","20 rue leblanc Paris","2 rue ambroise pare paris","1 place parvis notre dame paris",
              "149 rue de sevres paris","1 avenue claude vellafaux paris","184 faubourg saint antoine paris",
              "46 rue henri huchard paris")
hospital <- lapply(adresses,geocode)

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=47%20boulevard%20de%20l'hopital%20Paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=20%20rue%20leblanc%20Paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=2%20rue%20ambroise%20pare%20paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=1%20place%20parvis%20notre%20dame%20paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=149%20rue%20de%20sevres%20paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=1%20avenue%20claude%20vellafaux%20paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=184%20faubourg%20saint%20antoine%20paris&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=46%20rue%20henri%20huchard%20paris&sensor=false

and consider someone sick

library(ggmap)
adresses <- c("place Etoile paris")
sick <- lapply(adresses,geocode)

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=place%20Etoile%20paris&sensor=false

We can plot them on the map

darkcols=c("blue","red",rgb(1,0,0,.5))
plot(nodes[1:302,c("V3","V4")],pch=19,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",col="grey",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"],col="grey")
points(unlist(hospital)[names(unlist(hospital))=="lon"],unlist(hospital)[names(unlist(hospital))=="lat"],
       pch=forme[2],col=darkcols[2],cex=2)
points(unlist(sick)[names(unlist(sick))=="lon"],unlist(sick)[names(unlist(sick))=="lat"],
       pch=forme[1],col=darkcols[1],cex=2)

Let us move all those locations to the closest subsway station

close_knot=function(coord,knots=data.frame(X=nodes$V3,Y=nodes$V4,NO=1:nrow(nodes))){
  knots[which.min(distHaversine( coord[,c("lon","lat")] ,  knots[,c("X","Y")] , r =6378.137)),]
}
newhospital=data.frame(X=rep(NA,length(hospital)),Y=rep(NA,length(hospital)),NO=rep(NA,length(hospital)))
for(i in 1:length(hospital)) newhospital[i,]=as.numeric(close_knot(hospital[[i]]))
newsick=data.frame(X=rep(NA,length(sick)),Y=rep(NA,length(sick)),NO=rep(NA,length(sick)))
for(i in 1:length(sick)) newsick[i,]=as.numeric(close_knot(sick[[i]]))

plot(nodes[1:302,c("V3","V4")],pch=19,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",col="grey",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
                                nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"],col="grey")
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)

Using , it is possible to compute distance, for instance to get to the first hospital,

library(igraph)

## 
## Attaching package: 'igraph'

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

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

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

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

g <- make_graph(edges = as.vector(rbind(arcs[,1],arcs[,2])),directed=FALSE)
sp <- shortest_paths(g, newsick[1,"NO"], newhospital[1,"NO"], weights= arcs[,3])
sp

## $vpath
## $vpath[[1]]
## + 19/366 vertices, from 3c921b5:
##  [1] 318   7   8   9  10  11  12  13  14  15  87  86 362 242 243 244
## [17] 153 110 111
## 
## 
## $epath
## NULL
## 
## $predecessors
## NULL
## 
## $inbound_edges
## NULL

pth <- as.numeric(sp$vpath[[1]])

plot(nodes[,c("V3","V4")],pch=19,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",col="grey",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],
nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(nodes[pth,c("V3","V4")],col=darkcols[3],pch=19)
for(i in 1:(length(pth)-1)) segments(nodes[pth[i],"V3"],nodes[pth[i],"V4"],
nodes[pth[i+1],"V3"],nodes[pth[i+1],"V4"],col=darkcols[3],lwd=5)
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)

or the third one

library(igraph)
g <- make_graph(edges = as.vector(rbind(arcs[,1],arcs[,2])),directed=FALSE)
sp <- shortest_paths(g, newsick[1,"NO"], newhospital[3,"NO"], weights= arcs[,3])
pth <- as.numeric(sp$vpath[[1]])
plot(nodes[,c("V3","V4")],pch=19,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",col="grey",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(nodes[pth,c("V3","V4")],col=darkcols[3],pch=19)
for(i in 1:(length(pth)-1)) segments(nodes[pth[i],"V3"],nodes[pth[i],"V4"],
nodes[pth[i+1],"V3"],nodes[pth[i+1],"V4"],col=darkcols[3],lwd=5)
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)

Distances are respectively

distances(g, newsick[1,"NO"], newhospital[1,"NO"], weights= arcs[,3])

##      [,1]
## [1,] 7476

distances(g, newsick[1,"NO"], newhospital[3,"NO"], weights= arcs[,3])

##      [,1]
## [1,] 4438

On a graph, we get

plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE)
text(.125,.5,"DEMAND",srt=90,col=darkcols[1])
text(.875,.5,"SUPPLY",srt=270,col=darkcols[2])
agents=data.frame(X=c(newhospital$X,newsick$X),Y=c(newhospital$Y,newsick$Y),
                  Type=c(rep("Supply",nrow(newhospital)),c(rep("Demand",nrow(newsick)))))
i=which(agents$Type=="Demand")
y1=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
i=which(agents$Type=="Supply")
y2=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
id=which(agents$Type=="Demand")
for(i in 1:sum(agents$Type=="Demand")){
  d=rep(NA,sum(agents$Type=="Supply"))
  for(j in 1:sum(agents$Type=="Supply")){
    d[j]=distances(g, newsick[i,"NO"], newhospital[j,"NO"], weights= arcs[,3])
  }
  segments(.25,y1[i],.75,y2,col=c("black",darkcols[3])[1+(d==min(d))],lwd=c(1,6)[1+(d==min(d))])
  text(.5,(y2+y1[i])/2,round(d),pos=1,col=c("black",darkcols[3])[1+(d==min(d))])
}
i=which(agents$Type=="Demand")
points(rep(.25,length(i)),y1,
       pch=forme[1],col=darkcols[1],cex=3.2)
text(rep(.25,length(i)),y1,i,col="white")
i=which(agents$Type=="Supply")
points(rep(.75,length(i)),y2,
       pch=forme[2],col=darkcols[2],cex=3.2)
text(rep(.75,length(i)),y2,i,col="white")

Here is the closest hopital, for our sick agent

library(igraph)
g <- make_graph(edges = as.vector(rbind(arcs[,1],arcs[,2])),directed=FALSE)
sp <- shortest_paths(g, newsick[1,"NO"], newhospital[which.min(d),"NO"], weights= arcs[,3])
pth <- as.numeric(sp$vpath[[1]])

plot(nodes[,c("V3","V4")],pch=19,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",col="grey",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],
                                nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(nodes[pth,c("V3","V4")],col=darkcols[3],pch=19)
for(i in 1:(length(pth)-1)) segments(nodes[pth[i],"V3"],nodes[pth[i],"V4"],
                                nodes[pth[i+1],"V3"],nodes[pth[i+1],"V4"],col=darkcols[3],lwd=5)
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)

Voronoi sets

plot(nodes[,c("V3","V4")],pch=19,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",col="grey",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],
                                nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)
library(tripack)

## 
## Attaching package: 'tripack'

## The following objects are masked from 'package:igraph':
## 
##     convex.hull, triangles

loc=agents[agents$Type=="Supply",c("X","Y")]
names(loc)=c("x","y")
m=voronoi.mosaic(loc)
plot(m,add=TRUE)

affect1=function(x,y){
  which.min(distHaversine( c(x,y) ,  loc[,c("x","y")] , r =6378.137))
}
vx=seq(2.26,2.45,length=251)
vy=seq(48.8,48.9,length=251)
vz1=matrix(NA,251,251)
for(i in 1:251){
  for(j in 1:251){
    vz1[i,j]=affect1(vx[i],vy[j])
  } 
}
lightcols <- brewer.pal(8, "Pastel2")
image(vx,vy,vz1,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",axes=FALSE,col=lightcols)
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)

Now, if we use subway-based distances

close_knot=function(coord,knots=data.frame(X=nodes$V3,Y=nodes$V4,NO=1:nrow(nodes))){
  knots[which.min(distHaversine( coord[,c("lon","lat")] ,  knots[,c("X","Y")] , r =6378.137)),]
}

affect2=function(x,y){
pt=as.numeric(close_knot(data.frame(lon=x,lat=y)))
  for(j in 1:sum(agents$Type=="Supply")){
    d[j]=distances(g, pt[3], newhospital[j,"NO"], weights= arcs[,3])
  }
which.min(d)
}
vz2=matrix(NA,251,251)
for(i in 1:251){
  for(j in 1:251){
    vz2[i,j]=affect2(vx[i],vy[j])
  } 
}
image(vx,vy,vz2,xlim=c(2.26,2.45),ylim=c(48.8,48.9),xlab="",ylab="",axes=FALSE,col=lightcols)
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)