GoogleMaps + Bubble plot en R

En esta entrada muestro cómo graficar referencias geográficas empleando GoogleMaps y R. Para empezar necesitamos referencias geográficas que deseemos plotear, en mi caso graficaré unidades económicas de la Comisión Federal de Electricidad de los estados de Tlaxcala y Puebla, los datos los obtuve del Directorio Estadístico Nacional de Unidades Económicas (DENUE) del INEGI.

Instalamos y cargamos la librería

install.packages('RgoogleMaps')
library(RgoogleMaps)

Posterior a ello, definimos algunos parámetros como el zoom de nuestro mapa, la latitud y la longitud centro, y adicional a ello una distancia (Euclidiana) que será el radio de círculos concéntricos a partir de nuestro centro.

z0 <- 8                     # zoom
y0 <- 19.040442             # Latitud zero
x0 <- -98.250160            # longitud zero
k  <- 0.4                   # Define el salto del circulo
                            # Saltos de 11 Km por cada décima

Colocaremos en el centro del mapa (o cualquier coordenada), un marcador que nos indique alguna referencia adicional.

markers = paste0("&markers=color:blue|label:I|19.040442,-98.250160")

Acto seguido, obtener nuestro mapa base:

PuebMap  <- GetMap(center= c(y0,x0), zoom=z0, GRAYSCALE = TRUE,
             size=c(640,640),markers=markers,
             destfile = file.path(tempdir(),"interm.png"),
             maptype = 'mapmaker-hybrid')

Enseguida, importamos nuestra base de datos de las unidades económicas, la cual deberá tener una estructura parecida a esto:

> head(DF)
       id      lat       lon        x1 x2
1 3535207 18.87717 -97.72655   87.5724  3
2 3398690 18.49201 -97.40435 1167.6320 40
3 3263793 19.04774 -98.04557   87.5724  3
4 3467483 19.12102 -98.24399 1167.6320 40
5 3316191 18.95277 -97.91986   87.5724  3
6 3457112 19.06409 -98.23185   87.5724  3

x1 y x2 representan variables cuantitativas a representar mediante un gráfico de burbujas. Lo siguiente es elaborar el gráfico sobre nuestro mapa base.

tpm <- bubbleMap(DF,coords=c("lon","lat"),

                 crs=sp::CRS("+proj=longlat +datum=WGS84"),

                 max.radius=500,legendLoc="bottomright",

                 filename = file.path(tempdir(),"interm.png"),

                 map=PuebMap,zcol='zinc',

                 key.entries = c(100,250,500,750,2500,6000));

Lucirá de la siguiente manera:

Lo siguiente que haremos será graficar n círculos concentricos de radio n*k en nuestro mapa, con las siguientes líneas de comando:

r1 <- k*1

r2 <- k*2

r3 <- k*3

n  <- seq(0,359,0.05)

a  <- seq(1,360,0.05)

theta <- (360*n)/a

     x1 <- x0 + r1*cos(theta)

     y1 <- y0 + r1*sin(theta)

     x2 <- x0 + r2*cos(theta)

     y2 <- y0 + r2*sin(theta)

     x3 <- x0 + r3*cos(theta)

     y3 <- y0 + r3*sin(theta)

circle <- data.frame(c(x1,x2,x3),c(y1,y2,y3))

N <-length(x1)                        # Tamaño del vector

color <- c()                          # Inicializo vector

color[1:N]           <- 'gold'

color[(N+1):(2*N)]   <- 'darkorange'

color[(2*N+1):(3*N)] <- 'red'

{PlotOnStaticMap(PuebMap, lat = circle[,2], 

                lon = circle[,1], cex=0.5, pch=19,

                destfile = "interm.png", 

                col= color ,add=TRUE)};

Quedará algo como esto:

Podemos definir un nuevo nivel de zoom:

Por último y con fines de análisis geoespacial, para definir en que radio se encuentra cada una de las Unidades Económicas (UE), empleamos las siguientes líneas:

D <- dim(DF)[1]; D
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 1 o no
DF$PotR1 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r1^2){
      DF[j,'PotR1'] <- 'dentro'
      } else DF[j,'PotR1'] <- 'fuera'
}
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 2 o no
DF$PotR2 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r2^2){
      DF[j,'PotR2'] <- 'dentro'
      } else DF[j,'PotR2'] <- 'fuera'
}
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 3 o no
DF$PotR3 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r3^2){
      DF[j,'PotR3'] <- 'dentro'
      } else DF[j,'PotR3'] <- 'fuera'
}

Finalmente, en este ejemplo, se analiza a las variables cuantitativas en función de su pertenencia o no al área definida por cada radio, tal como el número de posibles empleados en cada radio:

sumER_1 <- 0
for (j in 1:D){
if (DF[j,'PotR1']=='dentro'){
    sumER_1 <- DF[j,'pers']+ sumER_1
    }
}
sumER_1                       # Potencial Radio 1

sumER_2 <- 0
for (j in 1:D){
if (DF[j,'PotR2']=='dentro'){
    sumER_2 <- DF[j,'pers']+ sumER_2
    }
}
sumER_2 <- sumER_2 - sumER_1             
sumER_2                       # Potencial Radio 2

sumER_3 <- 0
for (j in 1:D){
if (DF[j,'PotR3']=='dentro'){
    sumER_3 <- DF[j,'pers']+ sumER_3
    }
}
sumER_3 <- sumER_3 - sumER_2 - sumER_1   
sumER_3                       # Potencial Radio 3


sumER_x <- sum(DF['pers']) - sumER_1 - sumER_2 - sumER_3
sumER_x                       # Potencial fuera de cualquier radio

Información que finalmente representamos en una tabla:

> potencial <- data.frame(rbind(c(sumER_1, sumER_2, sumER_3, sumER_x)

+                               )

+                         )

> colnames(potencial) <- c(' Radio 1', ' Radio 2', ' Radio 3', ' Radio 4')

> rownames(potencial) <- 'Empleados'

> round(potencial)

           Radio 1  Radio 2  Radio 3  Radio 4

Empleados     1539      372      398      128

Lo que nos indica que existen aproximadamente 1539 empleados en el radio 1, que enfatizemos es de 44 Km a partir de la coordenada central que hemos fijado.