La prostitución y su entorno económico: un análisis espacial

Agrupación de puntos cercanos (polígonos) para mejorar la representación geoespacial

OSM

Prostitución

Clustering espacial

Centroide

Autor/a

Afiliación

Elena Vicente Carrasco

Universidad de Valencia

Fecha de publicación

30 de marzo de 2025

Input

Usando OpenStreetMap (OSM) se han estudiado y extraido datos sobre los establecimientos de cambio de divisa, burdeles, love hotels y licorerías con el fin de dar visibilidad a un tema marginal como la prostitución y los negocios que la rodean. Teniendo en cuenta la escasa información libre que podemos encontrar junto a la calidad de los puntos que obtenemos de OSM ha sido una labor difícil.

Dado que esta herramienta, OpenStreetMap, tiene un carácter colaborativo y de uso libre bajo una licencia abierta, hace que en muchas ocasiones no haya una estructura clara y estandarizada de los datos. De ahí que se hayan encontrado problemas a la hora de graficar los puntos. Así pues, se usarán técnicas para mejorar la calidad de los datos, no solo de las variables estudiadas, sino aplicable a cualquier conjunto de puntos.

Descripción

Se han detectado puntos que están demasiado cerca los unos de los otros como para ser edificios en sí, sino que representan el polígono que conforman esos edificios. Así pues se ha establecido una distancia limitante en metros para identificar si son edificios yuxtapuestos o simplemente puntos que delimitan el polígono del edificio. Para verlo de forma gráfica, cargaremos los datos obtenidos en prácticas anteriores sobre el territorio provincial de Alicante sobre la variable burdeles.

Cargamos las librerías necesarias:

library(tidyverse)
library(sf)
library(mapSpain)
library(leaflet)
library(shiny)
library(dbscan)
library(geosphere)

Leemos el conjunto de datos original generado en la Práctica 1 descargable en formato gpkg aqui. De este solo escogeremmos la información sobre burdeles como ya se ha comentado previamente:

antiguo_alc <- st_read('data/000141/alc.gpkg', layer ='Burdel')

Reading layer `Burdel' from data source 
  `C:\Users\scorp\OneDrive - Universitat de València\01 PROFESOR\06 Innovación Educativa\2023 2024\03 Web\datasets\data\000141\alc.gpkg' 
  using driver `GPKG'
Simple feature collection with 390 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.7369956 ymin: 37.99637 xmax: -0.0192765 ymax: 38.86795
Geodetic CRS:  WGS 84

nrow(antiguo_alc)

[1] 390

Este contiene 312 puntos que se pueden ver graficados a continuación:

antiguo_alc %>%
  leaflet::leaflet() %>%
  leaflet::addTiles()%>%  
  leaflet::addCircles(color = 'red', radius = 5)

Pero si nos acercamos a los puntos, vemos que ocurren cosas como estas:

Tratamiento

Así pues, se ha establecido una distancia mínima: 50 metros para agrupar los puntos que estén más cerca de esos 50 metros los unos de los otros y calcular un centroide. Es importante a tener en cuenta que para poder aplicar ese límite debemos extraer los datos de OpenStreetMap con crs 25830 cuyas unidades geográficas son en metros. Después se clusterizarán los puntos para así comprobar cuales son puntos únicos y cuales están demasiado cerca y con estos últimos calcular el punto medio.

df <- antiguo_alc
df$geom <- st_transform(df$geom, crs = 25830)

distancia <- 50 #distancia en metros
db_res <- dbscan(st_coordinates(df), eps = distancia, minPts = 2)
df$cluster <- db_res$cluster

knitr::kable(head(df, 20))

nombre	geom	cluster
Burdel	POINT (723106.5 4249929)	1
Burdel	POINT (723097.7 4249948)	1
Burdel	POINT (723102.7 4249950)	1
Burdel	POINT (723103.1 4249950)	1
Burdel	POINT (723103.6 4249950)	1
Burdel	POINT (723106.6 4249943)	1
Burdel	POINT (723104.5 4249942)	1
Burdel	POINT (723105.8 4249940)	1
Burdel	POINT (723113.8 4249943)	1
Burdel	POINT (723117.9 4249935)	1
Burdel	POINT (723105.6 4249931)	1
Burdel	POINT (723100.2 4249943)	1
Burdel	POINT (705198.3 4220584)	2
Burdel	POINT (705186 4220617)	2
Burdel	POINT (705207 4220611)	2
Burdel	POINT (705177 4220590)	2
Burdel	POINT (728026.6 4257802)	3
Burdel	POINT (728046.5 4257822)	3
Burdel	POINT (728054.7 4257814)	3
Burdel	POINT (728034.9 4257794)	3

duplicados <- df %>%
  filter(cluster != 0) %>%
  group_by(cluster)

knitr::kable(t(duplicados %>% distinct(duplicados$cluster) %>% select(cluster)))

cluster	1	2	3	4	5	6	7	8	9	10

nrow(duplicados)

[1] 390

Tal y como podemos comprobar en la tabla los puntos se agrupan en clusters del 1 al 10, en el caso de esta variable; y obtenemos un total de 312 duplicados. Después agrupamos y calculamos el centroide de los mismos para posteriormente juntarlos con los puntos únicos (en este caso ninguno). Por último, cambiamos el crs a 4326 para poder graficar y guardamos los datos.

punto_medio <- duplicados%>%
  summarise(nombre = first(nombre), geometry = st_centroid(st_union(geom))) %>%
  ungroup()%>%
  select(-cluster)
nrow(punto_medio)

[1] 10

unicos <- df %>%
  filter(cluster == 0) %>%
  select(-cluster)
  
union <- rbind(unicos, punto_medio)
nrow(union)

[1] 10

alc_4326 <- sf::st_transform(union, crs=4326)

suppressMessages(st_write(alc_4326, 'alc_burdel.gpkg', append = TRUE))

Updating layer `alc_burdel' to data source `alc_burdel.gpkg' using driver `GPKG'
Updating existing layer alc_burdel
Writing 10 features with 1 fields and geometry type Point.

En definitiva de los iniciales 312 “puntos” se han quedado en 10, algo mucho más razonable teniendo en cuenta los puntos que aparecen graficados en el mapa.

Output

Gracias a esto se ha podido obtener un mapa mucho más limpio y fiel a la realidad a pesar de que siguen apareciendo pocos establecimientos ya que muchos de estos no se registran como burdeles.

alc_4326 %>%
  leaflet::leaflet() %>%
  leaflet::addTiles()%>%  
  leaflet::addCircles(color = 'red', radius = 5, fillOpacity = 0.9)

Esto se puede aplicar a cualquier nivel geográfico y para las variables que decidamos. De hecho, si lo hacemos para todas las variables analizadas en una primera instancia tanto a nivel de España como de la Comunidad Valenciana o Alicante obtenemos los siguientes mapas:

[1] TRUE

Hacer clic para ver el código

# BOUNDING BOX
spain <- mapSpain::esp_get_country(moveCAN = FALSE)
spain <- st_transform(spain, crs = 4326)
bbox_spain <- st_bbox(c(xmin = -10, ymin = 35, xmax = 5, ymax = 45))

c_val <- mapSpain::esp_get_ccaa(ccaa = "Comunidad Valenciana")
c_val <- st_transform(c_val, crs = 4326)
bbox_c_val <- sf::st_bbox(c_val)

alc <- mapSpain::esp_get_prov(prov = "03")
alc <- st_transform(alc, crs = 4326)
bbox_alc <- sf::st_bbox(alc)

variables <- c("bureau_de_change", "brothel", "love_hotel", "alcohol")
label <- c("Cambio de divisa", "Burdel", "Hotel +18", "Licorería")

# OBTENER PUNTOS
get_osm_data <- function(bbox, loc, key, value, label) {
  q <- bbox %>%
    osmdata::opq(timeout = 10000) %>%
    osmdata::add_osm_feature(key = key, value = value)
  
  data <- osmdata::osmdata_sf(q)
  bbox_sfc <- st_as_sfc(bbox)
  points <- st_intersection(data$osm_points, loc) %>%
      mutate(nombre = label) %>%
      select(nombre, geometry)
  points <- st_transform(points, crs = 25830)
}

grafico <- function(bbox, loc) {
  name <- deparse(substitute(loc))
  # EXTRAEMOS LOS PUNTOS HACIENDO SO DE LA FUNCIÓN ANTERIOR
  exchange_points <- get_osm_data(bbox, loc, "amenity", variables[1], label[1])
  burdel_points <- get_osm_data(bbox, loc, "amenity", variables[2], label[2])
  love_hotel_points <- get_osm_data(bbox, loc, "amenity", variables[3], label[3])
  alcohol_points <- get_osm_data(bbox, loc, "shop", variables[4], label[4])
  
  df <- rbind(exchange_points, burdel_points, love_hotel_points, alcohol_points)
  
  distancia <- 50 #distancia en metros
  db_res <- dbscan(st_coordinates(df), eps = distancia, minPts = 2)
  df$cluster <- db_res$cluster
  
  duplicados <- df %>%
    filter(cluster != 0) %>%
    group_by(cluster)

  punto_medio <- duplicados%>%
    summarise(nombre = first(nombre), geometry = st_centroid(st_union(geometry))) %>%
    ungroup()%>%
    select(-cluster)
  
  unicos <- df %>%
    filter(cluster == 0) %>%
    select(-cluster)

  union <- rbind(unicos, punto_medio)
  
  df_4326 <- sf::st_transform(union, crs=4326)
  
  file_gpkg <- paste0(name, "_unico.gpkg")
  suppressWarnings(st_write(df_4326, file_gpkg, append = FALSE))
  
  # ESTABLECEMOS LOS COLORES PARA CADA CATEGORÍA
  colores <- c("red", "green", "orange", "purple")
  
  categorias <- unique(df_4326$nombre)
  pal <- colorFactor(colores[match(categorias, label)], 
                     domain = categorias)
  # MAPA BASE
  map <- leaflet() %>%
    addTiles()
  
  for (cat in categorias) { # esto es pq mfc solo tiene 1 variable 
    lista_variables <- df_4326 %>% filter(nombre == cat)
    
    if (nrow(lista_variables) > 0) { # MIENTRAS QUEDEN VARIABLES, AÑADIR PUNTOS
      map <- map %>%
        addCircles(
          data = lista_variables,
          color = pal(cat),
          fillOpacity = 1,
          label = ~nombre,
          popup = ~nombre,
          group = cat)
    }
  }
  
  # LEYENDA
  map <- map %>%
    addLegend(
      position = "bottomright",
      pal = pal,
      values = categorias,
      title = "Leyenda")
  # CAPAS
  map <- map %>%
    addLayersControl(
      overlayGroups = categorias,
      options = layersControlOptions(collapsed = T)) 
  # EXTRAS
  map %>%
    leaflet.extras::addResetMapButton() %>%
    leaflet.extras::addSearchOSM()
}

Mapa de España:

grafico(bbox_spain, spain)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Deleting layer `spain_unico' using driver `GPKG'
Writing layer `spain_unico' to data source `spain_unico.gpkg' using driver `GPKG'
Writing 811 features with 1 fields and geometry type Point.

Archivo generado con los puntos actualizados.

Mapa de la Comunidad Valenciana:

grafico(bbox_c_val, c_val)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Deleting layer `c_val_unico' using driver `GPKG'
Writing layer `c_val_unico' to data source `c_val_unico.gpkg' using driver `GPKG'
Writing 79 features with 1 fields and geometry type Point.

Archivo generado con los puntos actualizados.

Mapa de Alicante:

grafico(bbox_alc, alc)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Deleting layer `alc_unico' using driver `GPKG'
Writing layer `alc_unico' to data source `alc_unico.gpkg' using driver `GPKG'
Writing 48 features with 1 fields and geometry type Point.

Archivo generado con los puntos actualizados.

Proyectos de Innovación Educativa Emergente PIEE-2737007 y PIEE-3325394