Geoprocesamiento y Enriquecimiento de Datos Vitivinícolas: Análisis Urbanístico y Bioclimático en La Rioja

Autor/a
Afiliación

Adrián García Mirabet

Universidad de Valencia

Fecha de publicación

1 de abril de 2026

Input

Se ha obtenido un conjunto de datos espaciales correspondiente al Inventario de Elementos de los Barrios de Bodegas de La Rioja, extraído directamente de la Infraestructura de Datos Espaciales de La Rioja (IDERioja).

  • Fuente: Portal de datos abiertos del Gobierno de La Rioja (iderioja.larioja.org). Formato original: GeoJSON (estándar abierto de intercambio de datos geoespaciales). Naturaleza de los datos: El dataset original se distribuye mediante servicios de visualización de mapas (WFS/WMS) y descarga directa, proporcionando una capa de puntos que geolocaliza activos vitivinícolas tradicionales y modernos.

  • Estado inicial: El archivo presenta una estructura de metadatos orientada a la visualización web, incluyendo campos multimedia (URLs de imágenes y thumbnails) que, aunque valiosos para el usuario final, suponen una redundancia técnica para el análisis masivo de datos y el procesamiento estadístico.

  • Nota técnica: El dataset se distribuye originalmente en el sistema de referencia WGS 84 (EPSG:4326), lo que garantiza su compatibilidad con visores web globales, pero requiere una transformación para análisis métricos locales en España.

Descripción

Se detecta que el conjunto de datos original, aunque ofrece una base de geolocalización de altísimo valor por su escala micrométrica (geolocalizando elementos exactos en formato vectorial POINT), presenta severas limitaciones técnicas en bruto. Estas deficiencias impiden su uso directo para el análisis espacial avanzado y la investigación socioeconómica del “terroir” riojano:

  1. Ruido en la Tabla de Atributos: El archivo GeoJSON contiene un exceso de campos multimedia (URLs de fotografías y metadatos del visor web de IDERioja) que sobrecargan el fichero y la memoria sin aportar valor estadístico, haciendo indispensable un filtrado de variables.
  2. Fragmentación y Falta de Normalización: Existe una gran dispersión semántica en la variable “Tipo” (zarceras, tuferas, calados, bodegas…) sin una jerarquía patrimonial clara. Además, presenta inconsistencias tipográficas en los municipios que obligan a normalizar las cadenas de texto.
  3. Sistema de Referencia no Proyectado: El uso nativo de coordenadas geográficas WGS 84 (EPSG:4326) es el estándar web, pero imposibilita el cálculo métrico preciso de distancias topológicas y urbanísticas, requiriendo su transformación a ETRS89.
  4. Inexistencia de la Dimensión Vertical (Eje Z): El modelo es bidimensional. La carencia de datos de altitud impide estudiar la relación estructural entre el patrimonio excavado y la orografía del terreno.
  5. Carencia de Variables Bioclimáticas: El dataset aísla los activos de su entorno ambiental, careciendo de información térmica que permita relacionar la expansión del tejido industrial moderno con estrategias de adaptación frente al cambio climático.

Tratamiento

En primer lugar debemos cargar las librerías necesarias para ejecutar el código.

Para asegurarnos de que todas las librerías estén disponibles, utilizamos el siguiente código que las carga y, en caso de que falte alguna, la instala automáticamente.

packages <- c("sf", "dplyr", "elevatr", "leaflet", "geodata", "terra","RColorBrewer")

for (p in packages) {
  if (!requireNamespace(p, quietly = TRUE)) {
    # Añadimos el 'repos' para evitar el error de CRAN al renderizar en Quarto
    install.packages(p, dependencies = TRUE, repos = "http://cran.us.r-project.org")
  }
  library(p, character.only = TRUE)
}

Ahora cargaremos los datos oficiales de la rioja con el siguiente código:

datos_oficiales_geo <- st_read("../data/2526020019/Elementos de los Barrios de Bodegas.geojson", quiet=T)

Y ahora seleccionaremos las variables que nos interesan de este dataset limpiando las variables de tipo imagen etc que no son del todo interesantes para nuestro estudio:

datos_vinedo <- datos_oficiales_geo %>%
  select(id, Nombre, Tipo, Municipio, geometry) %>%
  mutate(Municipio = toupper(Municipio)) # Normalización para evitar duplicados por tildes/mayúsculas

Y ahora pasarermos a transformar el CRS de nuestro dataset a ETRS89 que es el oficial en España:

datos_vinedo <- datos_vinedo %>%
  st_transform(crs = 25830)

Comenzamos a enriquecer la base de datos donde vamos a añadir las columnas “subzona” y “categoria” para obtener las zonas en las que se encuentran las bodejas segun la zona (alta, oriental o media). En la parte de categoría es muy interesante obtener de la variable tipo si es patrimonio histórico (calados) o es de actividad industrial, que se correponde a bodega. De esta manera podremos comparar el patrimonio histórico con la actividad industrial.

datos_finales <- datos_vinedo %>%
  mutate(
    Subzona = case_when(
      Municipio %in% c("HARO", "BRIONES", "CENICERO", "FUENMAYOR") ~ "Rioja Alta",
      Municipio %in% c("ALFARO", "ALDEANUEVA DE EBRO", "CALAHORRA") ~ "Rioja Oriental",
      TRUE ~ "Otras zonas / Rioja Media"
    ),
    Categoria = case_when(
      Tipo == "Calado" ~ "Patrimonio Histórico",
      Tipo == "Bodega" ~ "Activo Industrial",
      TRUE ~ "Otros"
    ),
    UTM_X = st_coordinates(.)[,1],
    UTM_Y = st_coordinates(.)[,2]
  )




table(datos_finales$Tipo)

Acceso calado       Agrario      Aguadojo       Almacén        Bodega 
            7             6             4             1          1664 
       Calado          Caño      Castillo      Cubachón          Lago 
         3014           108             3             4             4 
        Otros       Palacio     Religioso   Residencial   Respiradero 
           45             1             7             2            19 
        Solar        Tufera  Vitivinícola        Yesera       Zarcera 
            5           406            46             1             9

Al observar los valores que toma la variable tipo redifinimos la variable segun los nuevos valores de imputación que determinamos creando los valores “Espacios de crianza”, “Elementos de ventilacion”, “Infraestructura industrial”, “Arquitectura Singular” de esta manera obtenemos una base de datos más refinada.

datos_finales <- datos_finales %>%
  mutate(
    Categoria_Patrimonial = case_when(
      Tipo %in% c("Calado", "Acceso calado", "Caño") ~ "Espacios de Crianza",
      Tipo %in% c("Tufera", "Respiradero", "Zarcera") ~ "Elementos de Ventilación",
      Tipo == "Bodega" ~ "Infraestructura Industrial",
      Tipo %in% c("Castillo", "Palacio", "Religioso") ~ "Arquitectura Singular",
      TRUE ~ "Otros Elementos Auxiliares"
    )
  )

Altitud

Nuestra mayor mejora para la base de datos es que mediante la librería “elevater” vamos a ser capaces de añadir una nueva variable para poder comprender la distribución de las bodegas según la altitud de las mismas, lo cual es muy interesante al hacer estudios de cambio climático y temperatura que vamos a poder observar con posterioridad. Usaremos un zoom de 12 con un buen detalle. De esta manera obtenemos la nueva variable “elevation”

datos_finales_z <- get_elev_point(datos_finales, 
                                  prj = st_crs(datos_finales), 
                                  src = "aws", 
                                  z = 12)
Mosaicing & Projecting
Warning: [extract] transforming vector data to the CRS of the raster
Note: Elevation units are in meters

Hemos hecho una pequeña comparación de la altitud media por categoría para observar las distribuciones de las bodegas o calados. De esta manera se puede determinar a que altitud necesita la uva ser cosechada.

analisis_altitud <- datos_finales_z %>%
  st_drop_geometry() %>%
  group_by(Categoria) %>%
  summarise(Altitud_Media = mean(elevation, na.rm = TRUE),
            Altitud_Min = min(elevation, na.rm = TRUE),
            Altitud_Max = max(elevation, na.rm = TRUE))

print(analisis_altitud)
# A tibble: 3 × 4
  Categoria            Altitud_Media Altitud_Min Altitud_Max
  <chr>                        <dbl>       <dbl>       <dbl>
1 Activo Industrial             546.         362         927
2 Otros                         552.         373        1109
3 Patrimonio Histórico          555.         329         928

De esta manera determinamos que al rededor de los 550 metros es donde se condensa la mayoría de la actividad vitivinícola en La Rioja.

Representación gráfica:

Mediante la librería leaflet vamos a realizar un estudio visual de la distribuciñin de las bodega según la altitud en la que se encientran para poder determinar las distribuciones de las mismas. Se necesita cambiar el crs para que leaflet funcione de manera correcta.

mapa_web <- datos_finales_z %>% st_transform(4326)

pal <- colorNumeric(palette = "YlOrRd", domain = mapa_web$elevation)

leaflet(mapa_web) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    color = ~pal(elevation),
    radius = 5,
    stroke = FALSE, 
    fillOpacity = 0.8,
    popup = ~paste0("<b>", Nombre, "</b><br>",
                    "Tipo: ", Tipo, "<br>",
                    "Altitud: ", round(elevation, 0), " m")
  ) %>%
  addLegend(pal = pal, values = ~elevation, title = "Altitud (m)", position = "bottomright")

Se observa una clara distribución que podría asimilarse a una linea con pendiente negativa que seguramente se corresponda a alguna característica del terreno que en estos momentos desconocemos pero en un futuro haremos un mejor análisis cartográfico.

ANÁLISIS DE PROXIMIDAD: TRADICIÓN VS INDUSTRIA

Vamos a tratar de entender mediante topología espacial avanzada tratar de entender el análisis cuantificando el éxodo urbanistico , vamos a tratar de entender si en estos momentos los calados y las bodegas siguen pegados como tradicionalmente o si se encuentran lejos, lo que vamos a estudiar es medir la distancia entre cada bodega y su calado más cercano. De esta manera podemos entender el modelo de negocio y clasificarlo numéricamente. Primero vamos a separar el dataset en patrimonio e industria

calados <- datos_finales_z %>% filter(Categoria == "Patrimonio Histórico")
bodegas <- datos_finales_z %>% filter(Categoria == "Activo Industrial")

Mediante la función st_distance calculamos la matriz de distancias de cada bodega a todos los calados. Posteriormente extraemos la sidtancia mínima de cada bodega a un calado. Mediante el case_when vamos a clasificar la distancia de cada bodega y calado según si se encuentra a menos de 150 metros, entre 150 metros y 1000 metros o si está muy separado. De esta manera podemos determinar lo cercano que está de cada una de los barrios históricos y calados

distancias_matriz <- st_distance(bodegas, calados)

bodegas <- bodegas %>%
  mutate(
    dist_al_patrimonio_m = as.numeric(apply(distancias_matriz, 1, min)),
    # Clasificamos el modelo de negocio según su distancia al barrio histórico
    Modelo_Urbanistico = case_when(
      dist_al_patrimonio_m < 150 ~ "Integrada en Barrio Histórico (<150m)",
      dist_al_patrimonio_m >= 150 & dist_al_patrimonio_m < 1000 ~ "Expansión Periférica (150m-1km)",
      TRUE ~ "Finca Aislada"
    )
  )

Ahora podemos realizar un resumen de la altitud media, el total de bodegas que hay según el tipo de modelo urbanístico en los que se encuentran las bodegas.

# Resumen para nuestras conclusiones
resumen_urbanistico <- bodegas %>%
  st_drop_geometry() %>%
  group_by(Modelo_Urbanistico) %>%
  summarise(Total_Bodegas = n(),
            Altitud_Media = mean(elevation, na.rm = TRUE)) %>%
  mutate(Porcentaje = round((Total_Bodegas / sum(Total_Bodegas)) * 100, 1))

print(resumen_urbanistico)
# A tibble: 3 × 4
  Modelo_Urbanistico                    Total_Bodegas Altitud_Media Porcentaje
  <chr>                                         <int>         <dbl>      <dbl>
1 Expansión Periférica (150m-1km)                 121          530.        7.3
2 Finca Aislada                                    26          707.        1.6
3 Integrada en Barrio Histórico (<150m)          1517          545.       91.2

Visualización Modelo urbanístico

Vamos a visualizar las modificaciones del dataset que hemos realizado con anterioridad, lo primero será cambiar el crs a wgs84 y a partir de ahí creamos una paleta según el tipo que se podrá ver en la leyenda.

Para crear el mapa interactivo lo dividimos en dos capas, la primera de calles, la segunda de vista satélite. Tambien dividimos en dos capas según si es calado histórico o bodega donde cada uno tendrá una forma y un color, donde el color de las bodegas estará dividido según el modelo urbanístico generado en el apartado anterior.

calados_wgs84 <- st_transform(calados, 4326)
bodegas_wgs84 <- st_transform(bodegas, 4326)

pal_urbana <- colorFactor(
  palette = c("#2ca25f", "#fdae61", "#d7191c"), 
  domain = bodegas_wgs84$Modelo_Urbanistico
)
mapa_proximidad_interactivo <- leaflet() %>%
  # Capa de fondo 1
  addProviderTiles(providers$CartoDB.Positron, group = "Mapa Callejero") %>%
  # Capa de fondo 2
  addProviderTiles(providers$Esri.WorldImagery, group = "Satélite") %>%
  
  # Capa 1
  addCircleMarkers(
    data = calados_wgs84,
    color = "grey",
    radius = 2,
    stroke = FALSE,
    fillOpacity = 0.5,
    group = "Calados Históricos"
  ) %>%
  
  # Capa 2
  addCircleMarkers(
    data = bodegas_wgs84,
    color = ~pal_urbana(Modelo_Urbanistico), # Color dinámico según el cálculo
    radius = 6,
    stroke = TRUE,
    weight = 1,
    fillOpacity = 0.9,
    group = "Bodegas Modernas",
    popup = ~paste0(
      "<div style='font-family: Arial, sans-serif; min-width: 200px;'>",
      "<b style='color:#333; font-size:14px;'>", Nombre, "</b><br>",
      "<hr style='margin: 5px 0;'>",
      "<b>Modelo:</b> <span style='color:", pal_urbana(Modelo_Urbanistico), ";'>", Modelo_Urbanistico, "</span><br>",
      "<b>Distancia al Patrimonio:</b> ", round(dist_al_patrimonio_m, 0), " m<br>",
      "<b>Altitud (Eje Z):</b> ", round(elevation, 0), " m",
      "</div>"
    )
  ) %>%
  
  addLegend(
    pal = pal_urbana,
    values = bodegas_wgs84$Modelo_Urbanistico,
    title = "Modelo Urbanístico",
    position = "bottomright"
  ) %>%
  
  addLayersControl(
    baseGroups = c("Mapa Callejero", "Satélite"),
    overlayGroups = c("Bodegas Modernas", "Calados Históricos"),
    options = layersControlOptions(collapsed = FALSE)
  )

#Visualización
mapa_proximidad_interactivo

Enriquecimiento bioclimático:Temperatura

Para continuar añadiendo valor e información a nuestros datasets vamos a obtener la temperatura media de cada bodega durante el año. Para realizar este proceso vamos a usar las librerias “geodata” y “terra” , vamos a preparar los puntos mediante la obtención de datos raster para España “ESP” y descargaremos las variables bioclimáticas “bio” , vamos a usar una resolución de 2.5 minutos que corrersponden a un arco de aproximadamente 4,5 kilometros. Y vamos a extraer los valores para los puntos donde nos quedamos con la capa BI01 que es la temperatura media anual, posteriormente con terra:extract asignamos a cada bodega la temperatura del pixel donde cae.

raster_climatico_spain <- worldclim_country(country = "ESP", var = "bio", res = 2.5, path = tempdir())
raster_temp_media <- raster_climatico_spain[[1]] 
names(raster_temp_media) <- "Temp_Media_Anual"
temp_extraida <- terra::extract(raster_temp_media, bodegas_wgs84)
bodegas_finalisimas <- bodegas_wgs84 %>%
  mutate(
    Temp_Media_C = temp_extraida$Temp_Media_Anual 
  )

Ahora procedemos a graficar las bodegas con esta nueva variable de temperatura media:

pal_temp <- colorNumeric(
  palette = "YlOrRd", 
  domain = bodegas_finalisimas$Temp_Media_C
)
mapa_temperatura <- leaflet(bodegas_finalisimas) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  
  addCircleMarkers(
    radius = 6,
    fillColor = ~pal_temp(Temp_Media_C), # Color según la temperatura
    color = "white",                     
    weight = 1,
    fillOpacity = 0.9,
    stroke = TRUE,
    
    popup = ~paste0(
      "<div style='font-family: sans-serif;'>",
      "<b>", Nombre, "</b><br>",
      "<hr style='margin:5px 0;'>",
      "<b>Temperatura Media:</b> <span style='color:red;'>", round(Temp_Media_C, 1), " ºC</span><br>",
      "<b>Altitud:</b> ", round(elevation, 0), " m<br>",
      "<b>Modelo:</b> ", Modelo_Urbanistico,
      "</div>"
    )
  ) %>%
  
  addLegend(
    pal = pal_temp, 
    values = ~Temp_Media_C, 
    title = "Temp. Media (ºC)", 
    position = "bottomright"
  )

mapa_temperatura

Ahora procedemos a realizar un informe final donde se puede observar según el modelo urbanístico que hemos definido en las primeras secciones donde vamos a observar el total, la altitud media de cada variable y la temperatura media de las mismas.

library(sf)
analisis_nicho_climatico <- bodegas_finalisimas %>%
  st_drop_geometry() %>%
  group_by(Modelo_Urbanistico) %>%
  summarise(
    Total = n(),
    Altitud_Media = mean(elevation, na.rm = TRUE),
    Temp_Media = mean(Temp_Media_C, na.rm = TRUE)
  ) %>%
  arrange(desc(Altitud_Media))

print(analisis_nicho_climatico)
# A tibble: 3 × 4
  Modelo_Urbanistico                    Total Altitud_Media Temp_Media
  <chr>                                 <int>         <dbl>      <dbl>
1 Finca Aislada                            26          707.       11.4
2 Integrada en Barrio Histórico (<150m)  1517          545.       12.2
3 Expansión Periférica (150m-1km)         121          530.       12.1

Output final integrado

Tras realizar todos los pasos anteriores integramos todo en un último dataset, donde observaremos en un gráfico final con toda la información integrada, integraremos el crs para comenzar, extraeremos los elementos auxiliares para que no se pierdan y el usuario pueda activarlos o no en el mapa

bodegas_mapa_final <- st_transform(bodegas_finalisimas, 4326)
calados_mapa_final <- st_transform(calados, 4326) # Contexto histórico


otros_elementos_mapa <- datos_finales_z %>% 
  filter(Categoria_Patrimonial %in% c("Elementos de Ventilación", "Arquitectura Singular", "Otros Elementos Auxiliares")) %>%
  st_transform(4326)

Ahora procedenos a realizar las paletas de color que luego se reflejarán en la leyenda . Usaremos el color gris para los calados históricos.

pal_urbana <- colorFactor(
  palette = c("#2ca25f", "#fdae61", "#d7191c"),
  domain = bodegas_mapa_final$Modelo_Urbanistico
)

Procedemos a crear el mapa final el cual va a aportar diferente información como la altitud de las bodegas que estamos estudiando, el tipo de bodega que es según la proximidad a calados, el municipio y la distancia a la que tiene patrimonio así como la temperatura media de la bodega en el año.

El código es largo y ligeramente complejo puesto a que se añade mucha información a este gráfico pero no deja de ser iterar lo mismo unas cuantas veces

mapa_resiliencia_bioclimatica <- leaflet() %>%
  # Capas de Fondo 
  addProviderTiles(providers$CartoDB.Positron, group = "Mapa Callejero") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "Satélite") %>%
  addProviderTiles(providers$OpenTopoMap, group = "Relieve Topográfico") %>%
  
  # CAPA 1 :Calados
  addCircleMarkers(
    data = calados_mapa_final,
    color = "#888888", # Gris suave
    radius = 2.5,
    stroke = FALSE,
    fillOpacity = 0.6,
    group = "Calados Históricos",
    popup = ~paste0("<b>", Nombre, "</b><br>Altitud: ", round(elevation, 0), " m")
  ) %>%
  
  # CAPA 2 Auxiliares
  addCircleMarkers(
    data = otros_elementos_mapa,
    color = "#444444", # Gris oscuro
    radius = 2,
    stroke = FALSE,
    fillOpacity = 0.8,
    group = "Elementos Auxiliares (Tuferas, etc.)",
    popup = ~paste0("<b>", Nombre, "</b><br>Tipo: ", Tipo)
  ) %>%
  
  # CAPA 3Bodegas Industriales y Análisis Bioclimático
  addCircleMarkers(
    data = bodegas_mapa_final,
    fillColor = ~pal_urbana(Modelo_Urbanistico), 
    color = "white",                             
    radius = 6.5,
    stroke = TRUE,
    weight = 1.2,
    fillOpacity = 0.9,
    group = "Bodegas Modernas (Análisis)",
    
    # POPUP 
    popup = ~paste0(
      "<div style='font-family: Arial, sans-serif; min-width: 250px;'>",
      "<b style='color:#333; font-size:15px;'>", Nombre, "</b><br>",
      "<hr style='margin: 8px 0;'>",
      "<b>Municipio:</b> ", Municipio, "<br>",
      "<b>Tipo:</b> ", Tipo, "<br>",
      "<hr style='margin: 8px 0;'>",
      "<b>Modelo Urbanístico:</b> <span style='color:", pal_urbana(Modelo_Urbanistico), ";'><b>", Modelo_Urbanistico, "</b></span><br>",
      "<b>Distancia al Patrimonio:</b> ", round(dist_al_patrimonio_m, 0), " m<br>",
      "<b>Altitud (Eje Z):</b> <span style='color:#007bff;'>", round(elevation, 0), " m</span><br>",
      "<b>Temperatura Media Anual:</b> <span style='color:#dc3545;'>", round(Temp_Media_C, 1), " ºC</span>",
      "</div>"
    )
  ) %>%
  
  # LEYENDA 
  addLegend(
    pal = pal_urbana,
    values = bodegas_mapa_final$Modelo_Urbanistico,
    title = "Modelo Urbanístico (Proximidad)",
    position = "bottomright"
  ) %>%
  
  # CONTROL DE CAPAS 
  addLayersControl(
    baseGroups = c("Mapa Callejero", "Satélite", "Relieve Topográfico"),
    overlayGroups = c("Bodegas Modernas (Análisis)", "Calados Históricos", "Elementos Auxiliares (Tuferas, etc.)"),
    options = layersControlOptions(collapsed = FALSE)
  )

mapa_resiliencia_bioclimatica

Output

Se ha obtenido un conjunto de datos espaciales multidimensional que trasciende la mera geolocalización bidimensional del archivo oficial. A través del geoprocesamiento, el dataset original de IDERioja ha sido enriquecido con tres capas analíticas fundamentales:

  1. Dimensión Altimétrica (Eje Z): Extracción de la cota sobre el nivel del mar mediante Modelos Digitales de Elevación (elevatr).
  2. Topología Urbanística: Cálculo de la distancia euclidiana de Nearest Neighbor entre la industria moderna y el patrimonio histórico para categorizar el modelo de implantación territorial.
  3. Nicho Bioclimático: Incorporación de la Temperatura Media Anual (BIO1) mediante la extracción de valores de la red global WorldClim.

Este enriquecimiento permite concluir empíricamente que la industria vitivinícola de expansión no solo se aleja del casco histórico, sino que coloniza cotas altimétricas superiores, encontrando refugios térmicos con temperaturas medias anuales casi 1 ºC más frías. Esta evidencia cartografía la adaptación real del sector frente al cambio climático.

#| message: false
#| warning: false

st_write(calados_mapa_final, 
         "bodegas_rioja_enriquecido.gpkg", 
         layer = "patrimonio_historico", 
         append = FALSE, 
         quiet = TRUE)

st_write(bodegas_mapa_final, 
         "bodegas_rioja_enriquecido.gpkg", 
         layer = "industria_bioclimatica", 
         append = TRUE, 
         quiet = TRUE)

st_write(otros_elementos_mapa, 
         "bodegas_rioja_enriquecido.gpkg", 
         layer = "elementos_auxiliares", 
         append = TRUE, 
         quiet = TRUE)

El/los fichero(s) generados con este procedimiento/técnica/metodología se puede(n) descargar de aquí.



Proyecto de Innovación Educativa Emergente (PIEE-3898312)