LinkedIn

Mapeo con R desde un archivo .shape

Voy a leer el archivo shape, el cual esta guardado con el sistema de coordenadas WGS-84 que utiliza google.

# Listo la lista de archivos dentro del directorio y tomo solo aquel que tiene la terminacion shape.
archivo <- gsub("\\.shp","",subset(list.files(),str_detect(list.files(), "\\.shp")))

# Leo el archivo shape obtenido anteriormente.
shape <- readOGR(dsn=".",layer=archivo)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\20389047504\Desktop\BACKUP\Lucio\Pedido Marco\Gobierno en Calle - Febrero 2020 - Escucha\Comuna 1", layer: "Retiro_"
## with 60 features
## It has 6 fields

Tenemos un archivo con 60 manzanas la cuales queremos dividir en grupos de a dos.

# Voy convertir el archivo shape en un data frame para facilitar el trabajo.
shape_df <- as.data.frame(shape)
# Ordeno segun su coordenada en y 
shape_df <- shape_df[order(shape_df$centro_y), ]
# Se organizan los grupos de acuerdo a la cercania entre las manzanas
i=1
x_min=1
y_min=1

## Cantidad de clusters en lo que queremos dividirlo
cant_clust = nrow(shape)/2 # La cantidad de manzanas que tenemos segun la cantidad de manzanas que queremos en cada grupo.
k=1
shape_df$cluster=""
shape_clusters=shape_df


sum_min=1
while(k<(cant_clust+1)){
  for(j in 2:nrow(shape_df)){
    
    y <- abs(shape_df$centro_y[i] - shape_df$centro_y[j])
    x <- abs(shape_df$centro_x[i] - shape_df$centro_x[j])
    suma <- x+y
    if(suma<sum_min){
      sum_min <- suma
      shape_limit <- data.frame(rbind(shape_df[i,], shape_df[j,]))
    }  
    
  }
  shape_clusters$cluster[shape_clusters$sm==shape_limit$sm[1]] <- k
  shape_clusters$cluster[shape_clusters$sm==shape_limit$sm[2]] <- k
  k=k+1
  shape_df <- shape_df[shape_df$sm != shape_limit$sm,]
  shape_df <- shape_df[shape_df$sm != shape_limit$sm,]
  sum_min=1
}
# Le coloco a cada manzana el grupo o cluster correspondiente
shape <- shape[order(shape$centro_y), ]
shape$cluster <- shape_clusters$cluster

id_cluster <- as.data.frame(table(shape$cluster)) 
###
# Vamos armando cada mapa.
for (i in 1:nrow(id_cluster)) { # Hago el for por el total de cluster, es decir un mapa por cluster
  
  CL  <- subset(shape, as.character(shape$cluster)==as.character(id_cluster[i,1]))#Hago un subset por el cluster que estamos tratando
  
#data frame de poligono para graficar en ggplot
  CLF <- fortify(CL)
  
  
#mapa zona de las "X" manzanas

  box <- make_bbox(CLF$long, CLF$lat, f = 1.8) # Arma el bounding box, es decir solo el area que contiene a las manzanas a graficar 
  zoom <- calc_zoom(box) # Calcula el zoom optimo a utilizar para sacar el mapa de google.
  centroid <- c((box["left"]+box["right"])/2,(box["top"]+box["bottom"])/2) # Calcula el centroide a partir del BBox
  m <-  ggmap(get_googlemap(center = centroid, zoom = zoom))
  
  #mapa grupo
  mapita <- m + 
    
    geom_polygon(data = CLF,
                 aes(x = long, y = lat, group = group),
                 alpha = 0.5, col="turquoise",fill="turquoise") + # Agrego los poligonos de cada manzana
    
    geom_text(data = CL@data, aes(x = as.numeric(as.character(CL@data$centro_x)), 
                                  y = as.numeric(as.character(CL@data$centro_y)), label= CL@data$sm), size=4.7) + # Agrego los numeros de manzanas
    ggtitle (paste("Mapa", row.names(id_cluster)[i] , sep=" ")) + # Pongo el titulo del mapa y lo numero segun el cluster que estoy tratando
    labs(x="Comuna 1") +
    
    theme (plot.title= element_text(size=19.5, face="bold", vjust=1, lineheight=0.8),
           axis.line=element_blank(), 
           axis.text.y=element_blank(), axis.ticks=element_blank(),
           axis.title.x=element_text(size=14), axis.text.x=element_blank(),
           axis.title.y=element_blank(),
           panel.background=element_blank(), panel.border=element_blank(),
           panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
           plot.background=element_blank())
  
  plot(mapita)
}