Ejercicio 2

En este ejercicio en lugar de hacer un “copia-pega” del ejercicio visto en clase quise ir un poco más allá e investigar el uso de google maps en R.

Para eso existe el paquete RgoogleMaps:

library(knitr)
library(rmarkdown)
library(PKI)

## Loading required package: base64enc

library(packrat)
library(ggplot2)
library(ggmap)

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

library(RgoogleMaps)

Así mismo, tenemos que conectatnos a la API de google:

register_google(key = "AIzaSyBh8VFmf7FaA6d8PIADsTi8vVrVGEbE4Xg")

# this sets your google map permanently
register_google(key = "AIzaSyBh8VFmf7FaA6d8PIADsTi8vVrVGEbE4Xg", write = TRUE)

## Replacing old key (AIzaSyBh8VFmf7FaA6d8PIADsTi8vVrVGEbE4Xg) with new key in C:\Users\zapic\OneDrive\Documentos/.Renviron

has_google_key()

## [1] TRUE

google_key()

## [1] "AIzaSyBh8VFmf7FaA6d8PIADsTi8vVrVGEbE4Xg"

has_google_client()

## [1] FALSE

has_google_signature()

## [1] FALSE

scrub_key("key=d_5iD")

## [1] "key=xxx"

scrub_key("key=d_5iD", "[your \\1]")

## [1] "key=[your key]"

scrub_key("signature=d_5iD")

## [1] "signature=xxx"

scrub_key("client=a_5sS&signature=d_5iD")

## [1] "client=xxx&signature=xxx"

Vamos a hacer el ejemplo desde un restaurante. Aquí podemos ver cómo se muestra en el centro del mapa de google sin problema:

#Forma 1
 localizacion = "Restaurante"  
#Forma 2
localizacion=c(lon=-8.710249 ,lat=42.238533)

restaurante <- get_map(location=localizacion, zoom= 18, source="google", 
                 maptype="roadmap", crop=TRUE)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=42.238533,-8.710249&zoom=18&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx

ggmap(restaurante)

Existen distintas formas de visualización. En Satélite:

restaurante2 = get_map(location=localizacion, zoom= 19, 
                 source="google", maptype="hybrid", crop=TRUE)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=42.238533,-8.710249&zoom=19&size=640x640&scale=2&maptype=hybrid&language=en-EN&key=xxx

ggmap(restaurante2)

o en acuarela. Aquí las vias de mayor afluencia de tráfico coinciden con los colores más rojos:

restaurante3 = get_map(location=localizacion, zoom= 15, 
                 source="stamen", maptype="watercolor", crop=TRUE)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=42.238533,-8.710249&zoom=15&size=640x640&scale=2&maptype=terrain&key=xxx

## Source : http://tile.stamen.com/watercolor/15/15589/12133.jpg

## Source : http://tile.stamen.com/watercolor/15/15590/12133.jpg

## Source : http://tile.stamen.com/watercolor/15/15591/12133.jpg

## Source : http://tile.stamen.com/watercolor/15/15592/12133.jpg

## Source : http://tile.stamen.com/watercolor/15/15589/12134.jpg

## Source : http://tile.stamen.com/watercolor/15/15590/12134.jpg

## Source : http://tile.stamen.com/watercolor/15/15591/12134.jpg

## Source : http://tile.stamen.com/watercolor/15/15592/12134.jpg

## Source : http://tile.stamen.com/watercolor/15/15589/12135.jpg

## Source : http://tile.stamen.com/watercolor/15/15590/12135.jpg

## Source : http://tile.stamen.com/watercolor/15/15591/12135.jpg

## Source : http://tile.stamen.com/watercolor/15/15592/12135.jpg

## Source : http://tile.stamen.com/watercolor/15/15589/12136.jpg

## Source : http://tile.stamen.com/watercolor/15/15590/12136.jpg

## Source : http://tile.stamen.com/watercolor/15/15591/12136.jpg

## Source : http://tile.stamen.com/watercolor/15/15592/12136.jpg

ggmap(restaurante3)

#GEOLOCALIZACION

Usamos las mismas librerías nombradas en clase:

# Librerías
.packages = c("TSP","maps","sp","maptools","geosphere")
library(geosphere)
library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

library(rgdal)

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/zapic/OneDrive/Documentos/R/win-library/3.6/rgdal/gdal
##  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: C:/Users/zapic/OneDrive/Documentos/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:packrat':
## 
##     init

library(maps)
library(mapdata)
library(ggmap)
library(marmap)

## Registered S3 methods overwritten by 'adehabitatMA':
##   method                       from
##   print.SpatialPixelsDataFrame sp  
##   print.SpatialPixels          sp

## 
## Attaching package: 'marmap'

## The following object is masked from 'package:raster':
## 
##     as.raster

## The following object is masked from 'package:grDevices':
## 
##     as.raster

library(lattice)
library(TSP)

Lo que vamos a estudiar es una ruta de reparto desde el restaurante a 6 clientes distintos:

rutas_vigo = read.csv2("rutas_vigo2.csv") 

Aquí hay un error de escrituta. Aunque en el excel se escriba “lat” o “Lat” la primera columna siempre añade esto “ï..” no he conseguido arreglarlo…

#Cargamos las coordenadas de las ciudades en coords
coordenadas_vigo = cbind(rutas_vigo$long, rutas_vigo$lat)

#cargamos en spdf un dataframe con los puntos espaciales de las coordenadas
spdf_vigo = SpatialPointsDataFrame(coordenadas_vigo, rutas_vigo)

Estuidamos las distancias entre todos los clientes:

# calcular distancia real entre los puntos utilizando el método distHaversine que mide la distancia más corta entre dos puntos. Este método asume una tierra esférica, ignorando los efectos elipsoidales.
distancia_real_vigo<-distm(coordenadas_vigo,fun = distHaversine)
distancia_real_vigo

##           [,1]      [,2]     [,3]     [,4]     [,5]     [,6]
## [1,]    0.0000  360.6743 2313.327 3645.508 5460.918 1578.831
## [2,]  360.6743    0.0000 2029.261 3446.507 5484.261 1662.164
## [3,] 2313.3266 2029.2608    0.000 1638.901 4612.548 2097.559
## [4,] 3645.5078 3446.5067 1638.901    0.000 3407.657 2712.226
## [5,] 5460.9185 5484.2613 4612.548 3407.657    0.000 3891.346
## [6,] 1578.8306 1662.1638 2097.559 2712.226 3891.346    0.000

Podemos ver que estas distancias tienen sentido y están correctas.

#Necesitamos crear una matriz de distancias de tipo TSP para que funcione solve_TSP
tsp_vigo<-TSP(distancia_real_vigo)
tsp_vigo

## object of class 'TSP' 
## 6 cities (distance 'unknown')

Generemos nuestro mapa de vigo usando ggmap con la fuente de google:

vigo=c(lon=-8.719431 ,lat=42.229714)

Vigo_basemap = ggmap::get_map(location=vigo, zoom= 14, source="google", maptype="roadmap", crop=TRUE)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=42.229714,-8.719431&zoom=14&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx

ggmap(Vigo_basemap)

mapa con nombres de lugares

# mapa con nombres de lugares

plot(as(spdf_vigo, "Spatial"), axes=TRUE,xlim=c(-15,5),ylim=c(36,44))

plot(Vigo_basemap)
# se puede reducir el tamaño de los puntos reduciendo el valor de cex
points(spdf_vigo, pch=20, cex=3, col="black")
# añadir etiquetas de nombre de ciudad, modificar pos para poner texto encima, debajo, derecha o izquierad del punto
text(spdf_vigo$long,spdf_vigo$lat,labels=spdf_vigo$name,pos=4)

`

#HAY QUE CORRER TODO JUNTO PARA QUE SALGAN LOS PUNTOS

#Librerías
.packages = c("TSP","maps","sp","maptools","geosphere")
library(geosphere)
library(maptools)
library(rgdal)
library(raster)
library(maps)
library(mapdata)
library(ggmap)
library(marmap)
library(lattice)

# Validación
.inst <- .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0)
  suppressPackageStartupMessages(install.packages(.packages[!.inst]))

# Carga
suppressPackageStartupMessages(library("TSP"))
suppressPackageStartupMessages(library("sp"))
suppressPackageStartupMessages(library("maps"))
suppressPackageStartupMessages(library("maptools"))
suppressPackageStartupMessages(library("geosphere"))

rutas_vigo <- read.csv2("rutas_vigo2.csv")
rutas_vigo

##        lat      long        name
## 1 42.23853 -8.710249 Restaurante
## 2 42.23874 -8.714616   Cliente_1
## 3 42.22769 -8.734188   Cliente_2
## 4 42.21330 -8.738445   Cliente_3
## 5 42.18950 -8.712457   Cliente_4
## 6 42.22437 -8.709138   Cliente_5

#Cargamos las coordenadas de las ciudades en coordenadas
coordenadas_vigo = cbind(rutas_vigo$long, rutas_vigo$lat)
#cargamos en spdf un dataframe con los puntos espaciales de las coordenadas
spdf_vigo = SpatialPointsDataFrame(coordenadas_vigo, rutas_vigo)


# calcular distancia real entre los puntos utilizando el método distHaversine que mide la distancia más corta entre dos puntos. Este método asume una tierra esférica, ignorando los efectos elipsoidales.
distancia_real_vigo<-distm(coordenadas_vigo,fun = distHaversine)
distancia_real_vigo

##           [,1]      [,2]     [,3]     [,4]     [,5]     [,6]
## [1,]    0.0000  360.6743 2313.327 3645.508 5460.918 1578.831
## [2,]  360.6743    0.0000 2029.261 3446.507 5484.261 1662.164
## [3,] 2313.3266 2029.2608    0.000 1638.901 4612.548 2097.559
## [4,] 3645.5078 3446.5067 1638.901    0.000 3407.657 2712.226
## [5,] 5460.9185 5484.2613 4612.548 3407.657    0.000 3891.346
## [6,] 1578.8306 1662.1638 2097.559 2712.226 3891.346    0.000

#Necesitamos crear una matriz de distancias de tipo TSP para que funcione solve_TSP
tsp_vigo<-TSP(distancia_real_vigo)
tsp_vigo

## object of class 'TSP' 
## 6 cities (distance 'unknown')

#Descargamos el mapa de España, con level=1 lo tendremos por comunidad autónoma
vigo=c(lon=-8.719431 ,lat=42.229714)
Vigo_basemap = ggmap::get_map(location=vigo, zoom= 14, source="google", maptype="roadmap", crop=TRUE)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=42.229714,-8.719431&zoom=14&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx

# mapa con nombres de los destinos#plot(as(spdf_vigo, "Spatial"), axes=TRUE,xlim=c(-50,50),ylim=c(50,50))
#ggmap(Vigo_basemap)
#points(spdf_vigo, pch=10, cex=2, col="red")
#text(spdf_vigo$long,spdf_vigo$lat,labels=spdf_vigo$name,pos=4)

Podemos ver que nos da un error… no pone los puntos tal y donde debería.

Calculo de ruta usando Nearest-Neighbour Algorithm (NNA)

Como hemos dicho antes; vamos a calcular la ruta óptima para un repartidor del restaurante que tiene que atender a 6 clientes.

#Calculamos el tour con el método por defecto y decimos que empiece por el Cliente 1:
tour_vigo <- solve_TSP(tsp_vigo,method = "nn", start=1)

tour_vigo1 <- SpatialLines(list(Lines(list(Line(coordenadas_vigo[c(tour_vigo, tour_vigo[1]),])), ID="1")))


#Generamos el mapa otra vez:
#plot(as(spdf_vigo, "Spatial"), axes=TRUE,xlim=c(-10,5),ylim=c(36,44))
#ggmap(Vigo_basemap)

# se puede reducir el tamaño de los puntos reduciendo el valor de cex
#points(spdf_vigo, pch=20, cex=3, col="black")

# añadimos las etiquetas a cada destino:
#text(spdf_vigo$long,spdf_vigo$lat,labels=spdf_vigo$name,pos=4)

#Sacamos el tour más óptimo:
#plot(tour_vigo1, add=TRUE, col = "red")
#points(coordenadas_vigo, pch=3, cex=0.4, col="black")

tail(tour_vigo)

## 1 2 6 3 4 5 
## 1 2 6 3 4 5

tour_vigo

## object of class 'TOUR' 
## result of method 'nn' for 6 cities
## tour length: 14627.87

La ruta de reparto que nos proporciona es correcta:

\[ \begin{align*} Restaurante\Longrightarrow Cliente 2\Longrightarrow 6\Longrightarrow 3\Longrightarrow 4\Longrightarrow 5\end{align*} \]

En su momento conseguí que dibujase una ruta sobre el mapa de Vigo, donde los puntos no coincidian encima del mapa PERO la relación de situación y distancias entre ellos tenía sentido; pero algo toqué para tratar de mejorarlo y lo perdí sin poder recuperarlo.

Creo que el problema está en la lectura del CSV.

Aquí guardo una imagen de la ruta que podría salir.

png("rutas_vigo.png",

    width=2400,height=1620,units="px",

    pointsize=24,bg="white",res=300)

tour_vigo <- solve_TSP(tsp_vigo,method = "nn", start=1)
tour_vigo1 <- SpatialLines(list(Lines(list(Line(coordenadas_vigo[c(tour_vigo, tour_vigo[1]),])), ID="1")))
plot(as(spdf_vigo, "Spatial"), axes=TRUE,xlim=c(-10,5),ylim=c(36,44))
ggmap(Vigo_basemap)
points(spdf_vigo, pch=20, cex=3, col="black")
text(spdf_vigo$long,spdf_vigo$lat,labels=spdf_vigo$name,pos=4)
plot(tour_vigo1, add=TRUE, col = "red")
points(coordenadas_vigo, pch=3, cex=0.4, col="black")
tail(tour_vigo)

## 1 2 6 3 4 5 
## 1 2 6 3 4 5

dev.off()

## png 
##   2