Mapas del tesoro

Mapas del tesoro
true

Una vuelta más sobre el generador de islas aleatorias

Hace unas semanas publiqué un artículo dedicado a generar islas aeatorias. Fue un tema que gustó mucho y también a mis hijos, por lo que estoy realizando una versión en Shiny que permita ejecutar en línea el código y estará disponible para todos aquí.

También he visto que se puede aprovechar y darle alguna vuelta más profesional al asunto, en concreto para aprender de SIG (Sistemas de Información Geográfica), así que, manos a la obra que vamos a hacer mapas con R, vamos a crear Mapas del Tesoro.

Mapa del tesoro

Mapa del tesoro

Simple features

Las islas que generaba en el artículo anterior no tienen base geográfica, simplemente se pintan en un lienzo como polígonos. Usaremos las funciones de ese articulo y las representaremos como objetos espaciales con la librería más avanzada de SIG en R , con sf.

El paquete sf (simple features) que traducido es algo asi como objetos simples es el más avanzado y en él se integran de forma nativa los objetos geográficos en R. Las librerías que se usaban hasta ahora para representar mapas en R se basaban en otras librerías principalmente en sp (spatial), que era fantástica, pero adolece de muchos problemas que con sf se completo.

En sf existen funciones para todo lo que suponga operaciones y matemática geográfica vectorial, es decir para el manejo de objetos simples en intersecciones, uniones, buffer, transformaciones afines, simplificación, corte, transformación de tipo….

El problema es que sf es una librería compleja, que requiere de práctica y tiempo para comprender y no es sencillo de comprender al principio. Pero para eso estamos aquí, para poner ejemplos y práctica que nos ayuden.

Objetivo

El objetivo de este taller será obtener un Mapa del Tesoro de esos típicos de libro de aventuras, con una gran isla principal y alguna más pequeña en la costa cercana.

El mapa, para ser completo, debe tener tambien los nombres de algunos accidentes geográficos o puntos singulares como cabos, playas y por supuesto el punto donde está enterrado el tesoro.

Rescatamos el código de crear_isla

En el articulo anterior islas aeatorias programamos estas 5 funciones:

  • puntomedio, que calcula el punto medio entre dos puntos dados.
  • pol_cero, que crea un polígono aletorio a partir de los valores de diámetro medio y número de lados.
  • div_pol, función que divide cada lado del polígono de entrada en dos y genera distorsiones caóticas aleatorias distanciando el punto medio sobre la perpendicular al lado.
  • div_pol_n función recursiva que hace n veces la función anterior div_pol
  • crea_isla función que llama a las anteriores y crea el polígono de una isla aleatoria y lo pinta.

El esquema que tenemos en mente es el siguiente:

Tenemos que hacer algún retoque de la función crea_isla(), para que en lugar de pintar el resultado, devuelva un polígono en formato data.frame, eso os lo dejamos de tarea.

Por cierto que revisando el código he descubierto un fallo en la función div_pol, pues se duplicaban los puntos de coordenadas en cada división del poligono, ahora ya está corregido, ¡cosas de la programación rápida!.

Esquema de programación

Os muestro un esquema final con las funciones para hacer el mapa del tesoro. He tenido que meter bastantes funciones que no pensé al inicio por los problemas que surgen en la programación y cosas inesperadas que hay que ir solucionando.

He añadido una librería para manejar fuentes de texto ampliadas, lo que me permitirá que el mapa use tipografías de apariencia antigua…. es extrafont.

  # para usar otras caligrafias
  library(extrafont)
  #font_import() solo una vez tras instalar
  loadfonts() # carga llas fuentes

Como dije, he usado solamente la librería sf, pues tenemos en ella todas las funciones necesarias. Los objetos espaciales se almacenan en sf como en listas, dentro de la parte de geometría sfc, y como datos en el sf. Haré un articulo especifico para desentrañar esto más adelante.

sf tiene también la ventaja de usar el código enlazado con pipe (de tubo), usando la librería tidyverse.

Las principales funciones que usaremos sobre geometrias sf en este taller son: st_cast, st_centroid,st_buffer, st_combine, st_union, st_difference…

Funciones mapa tesoro 1

Siguiendo el esquema, las funciones que programamos para hacer el mapa del tesoro son las siguientes:

  • pol_to_sf()–> que devuelve un objeto sf a partir de un data.frame del polígono.
  • validar_contorno()–> esta función es muy importante, y es el resultado de muchos dolores de cabeza, pues los poligonos aleatorios son muchas veces extraños y tienen lados que se cruzan, lo que genera un MULTIPOLYGON o directamente geometrías no válidas. Con esta función nos quedamos unicamente un un polígono simple en fomrato sf.
  • multipol_a_uno()–> que coge el poligono mas grande de un MULTIPOLYGON

El código hasta aquí, es el siguiente:

################################################
# Funciones para hacer un mapa del tesoro con R
# Autor: Fernando Villalba Bergado
# Fecha: 2018-2019
################################################
  require(extrafont)
  library(tidyverse)
  library(sf)

  # Funcion pol_to_sf:
  # para transformar el poligono en objeto o capa espacial sf
  # por defecto coge el CRS epsg 25831 
    pol_to_sf<-function(pol,epsg=25831){
      capa <- pol %>%
        st_as_sf(coords = c("x", "y"), crs = epsg) %>%
        summarise(geometry = st_combine(geometry)) %>%
        st_cast("POLYGON")
      st_crs(capa)<-epsg
      return(capa)
    }

  # Funcion validar_contorno: 
  # sirve para corregir fallos en los poligonos en los que cruzan lados
    validar_contorno<-function(pol){
      str1<-st_crs(pol)
      if(is.na(str1)){str1=25831}
  #    pol<-st_cast(pol,"POLYGON")
      if(!st_is_valid(pol)){
        pol<-st_buffer(pol, 10)# crea un buffer
        pol<-st_combine(pol)
        st_crs(pol)<-str1#25831
      }
      # comprueba que no es multipol
      if(class(pol)[[1]]=="sfc_MULTIPOLYGON"){
        pol<-multipol_a_uno(pol)
        st_crs(pol)<-str1
      }  
      return(pol)
    }

      
  ## Función que coje solo el poligono de mayor área de un multipol
    multipol_a_uno<-function(mpol){
       aux1<-st_cast(mpol,"POLYGON")
    #  aux1<-st_cast(aq,"POLYGON")
        areas<-st_area(aux1)
        aux2<-data.frame(area=areas,id=1:length(areas))
        #selecciona la de mayor area
        solo1<-head(aux2[order(aux2$area,decreasing = T),2],1)
        # selecciona solo la de mayor area
        aux3<-st_sf(st_geometry(aux1)[solo1])
        return(aux3)
    }

# prueba de las funciones
isla_prueba<-crea_isla(R=5000)%>% pol_to_sf()%>%validar_contorno()
plot(isla_prueba,axes=T,
     main="Polígono-isla generada aleatoriamente",
     graticule=T,bg="aliceblue",col="antiquewhite4",
     family="Calibri") 

Funciones mapa tesoro 2

Con las funciones anteriores generamos un polígono con forma de isla y lo convertimos en un objeto geográfico con caracteristicas SIG, pero para crear un mapa del tesoro tenemos que currar un poco mas.

Así que creamos más funciones siguendo el esquema inicial:

  • isla_chica() que genera una isla en el contorno de una isla mayor, para esto usaremos la función buffer de sf.
  • archipielago() que genera varias islas chicas, llamando a la funcion anterior varias veces.
  • junta_igran_archi() que une la isla grande y las islas del archipielago
  • nom_islas() que da nombre a las islas generadas
  • p_interes()que genera puntos aleatorios de caracteristicas concretas para la geografia de las islas
  • pinta_mapa_negro() función que pinta el mapa en tonos negros

El código detallado es el siguiente:

################################################
# Funciones para hacer un mapa del tesoro con R
# Autor: Fernando Villalba Bergado
# Fecha: 2018-2019
################################################
  require(extrafont)
  require(prettymapr) 
  require(tidyverse)
  require(sf)

  # 2. Función que crear una isla chica de radio r 
  #  y la posiciona en la costa de la isla grande
    isla_chica<-function(islagrande,r){
        ichi<-crea_isla(R=r)%>% pol_to_sf()%>%validar_contorno()
        # Posicion de la isla en el bufer
        posi_ichi <- st_buffer(st_geometry(islagrande),r*2)
        # Mueve la isla peq a un punto aleatorio del buffer
        n<-dim(st_coordinates(posi_ichi))  
        dado<-as.integer(runif(1,1,n))# num aleatorio entre la dim de dado
        # vector desplazamiento
        x_chica<-st_coordinates(posi_ichi)[dado,1]
        y_chica<-st_coordinates(posi_ichi)[dado,2]
        # nuevo punto donde mover la isla chica
        p = st_point(c(x_chica,y_chica))
        #convierto el punto en geometria sfc
        p<-st_sfc(p, crs = 25831)
        #convierto la geometria en sf
        p<-st_sf(p)
        #sumo las geometrias sf es decir, 
        # desplazamos el poligono al punto
        aux<-ichi+p
        # por si acaso convierto en sf
        aux<-st_sf(aux,crs = 25831)
        return(aux)
    }  
  

  # 3. funcion que generea varias islas, un archipielago
  # devuelve una lista con las islas y el numero de estas  
    archipielago<-function(zona,radio=5000,n_islas=NA){ 
    # generamos unas islas chicas 
     n_islas<-ifelse(is.na(n_islas),as.integer(runif(1,2,5)),n_islas )
     archi<-list()# vacia la lista
        # crea el poligono de buffer al que va sumando islas
        todas_las_islas<-zona # el primero con la isla grande
        #todas_las_islas<-validar_contorno(todas_las_islas)
        # bucle que crea las islas
        for(i in 1:n_islas){
          # calucla un radio aleatorio de la isla chica
          # como media 1/10 de la grande
          r_isla<-radio/(runif(1,8,11))
          # calcula el bufer del contorno
          contorno<-st_buffer(todas_las_islas,(r_isla+500))
          # Llama a la función isla_chica
          archi[[i]]<-isla_chica(contorno,r_isla) %>% st_sf()
          #valida el contorno
          archi[[i]]<-validar_contorno(archi[[i]])          
          #une las islas para crear el nuevo contorno
          todas_las_islas<-st_union(todas_las_islas,archi[[i]])
        }
        # devuelve una lista de 2
        #aux<-list(archi,n_islas)
        return(archi)
    }  
  

    # Funcion junta_igran_archi
    # junta la isla principal y el archipielago de islas chicas
     junta_igran_archi<-function(igrande, archi){
       #archi<-a
       #igrande<-islagrande
       n_islas<-length(archi)
          # añado a la lista la islagrade
            archi[[n_islas+1]]<-igrande
            # unimos todas las islas de la lista archi 
            islas <-purrr::reduce(archi,st_union)
            # Dividimos en poligonos cerrados
            islas<-st_cast(islas,"POLYGON")
      return(islas)
     }

   
   # Funcion nom_islas:
   # da nombre a las islas aleatorias de entre los de la lista de la funcion
   nom_islas<-function(n_islas){
      nombres_islas<-c("Cabeza Rota","I. Pájaros",
                       "Isla Chica", "Islote de los Naúfragos",
                       "I. Esperanza", "I. Capitán Vilber",
                       "Isla Soledad","Peñón de la Negra",
                       "Punta del Moro")
      # seleccion de los nombres anteriores unos aleatorios
      # según el número de islitas
      nom_I<-sample(nombres_islas,n_islas)
      # Añadimos el nombre de la isla grande
      nom_I<-c(nom_I,"Isla Perdida" )  
      return(nom_I)
   }  
   
    # Funcion  p_interes:
    # Crea una capa de puntos de interes en la isla.
       # 2 puntos vip
       # 1 playa
       # 2 cabos
    # la funcion devuelve una lista con 3 objetos:
    #   puntos_vip,puntos_cabo,ruta_pirata
    p_interes<-function(islagrande){
      # calculo la linea de costa isla grande
      lin_costa<-st_cast(islagrande,"LINESTRING")
      #borra los puntos
      puntos_vip<-NULL
      puntos_cabo<-NULL
      punto_1<-NULL
      # Puntos interiores 
          npuntos<-2
          # ojo que esta funcion st_sample muchas veces no devuelve el num especificado en size,
          # por eso voy a calcular 2 mas y me quedo con los 2 primeros.
          puntos_vip<-st_sample(st_buffer(islagrande,-200),size=npuntos+3)
          puntos_vip<-st_sf(puntos_vip)
          puntos_vip<-puntos_vip[1:npuntos,]
      # Punto enlinea de costa para playa
          npuntos<-1
          punto_1<-st_sample(lin_costa,size=npuntos+2,type="random")
          punto_1<-st_sf(punto_1)
          punto_1<-st_cast(punto_1,"POINT")
          punto_1<-punto_1[1,] # por quedarme solo con 1
      # Unimos los 3
          puntos_vip<-st_union(punto_1,puntos_vip)
          puntos_vip<-st_difference(st_cast(puntos_vip,"POINT"))
      # Asignamos los nombres de los puntos      
          puntos_vip$nombre<-c("Playa Mosquito","Palmera Alta","Tesoro")
          
      # Cabos
          # simplificamos el contorno mucho para identificar vertices
          cape<-st_simplify(islagrande,dTolerance =2000)#(radio-1000))
          # convierto a linea
          cape<-st_cast(cape,"LINESTRING")
          # convierto a puntos
          cape<-st_cast(cape,"POINT")
          #st_crs(cape)
          npuntos<-2
          puntos_cabo<-st_sample(cape,size=npuntos+2,type="random")
          puntos_cabo<-st_cast(puntos_cabo,"POINT")
          puntos_cabo<-st_sf(puntos_cabo)
          puntos_cabo<-puntos_cabo[1:npuntos,] # por quedarme solo con 2
          puntos_cabo$nombre<-sample(c("Cabo Tormentas", "Cabo de la Desesperanza",
                                       "Punta Negra","Punta de las flechas"),npuntos)
        # creo una linea on la ruta del tesoro  
          ruta_pirata<-st_cast(st_combine(puntos_vip), "LINESTRING")
          aux<-list(puntos_vip,puntos_cabo,ruta_pirata)
          return(aux)
    }   

    
    # Funcion pinta_mapa:
    #   pinta el mapa completo del tesoro en plot base
    pinta_mapa_negro<-function(is2,nom_I,p_int){
        #Calcula el centroide de cada isla
        is2_c<-st_centroid(is2,of_largest_polygon = FALSE)
        
        # Pintamos todo  
          # pinta las islas
            plot(st_geometry(is2), col="black", border = "black",  axes = TRUE, main="Mapa del tesoro del capitán Vilber", bg="gray81",family="Old English Text MT")  #darkslategray1
          # Pinta los centroides
            plot(st_geometry(st_centroid(is2)),pch = 3, col = 'white', add = TRUE)
          # Pinta las lineas e costa del mar
            plot(st_cast(st_buffer(is2,50),"MULTILINESTRING"), col = 'gray20', add = TRUE)
            plot(st_cast(st_buffer(is2,100),"MULTILINESTRING"), col = 'gray23', add = TRUE)
            plot(st_cast(st_buffer(is2,300),"MULTILINESTRING"), col = 'gray29', add = TRUE)
            plot(st_cast(st_buffer(is2,500),"MULTILINESTRING"), col = 'gray32', add = TRUE)
        # Etiquetas de las islas
        text(st_coordinates(is2_c), labels = nom_I, pos = 4,family="Old English Text MT",col = 'white')
        # Puntos de interes
          plot(p_int[[1]],add=T,pch=4,cex=2,lwd=2,col="gray42")
          plot(p_int[[2]],add=T,col="red",pch=4,cex=2)
          plot(p_int[[3]],add=T, col="red",lwd=2,lty=2)
        #Etiquetas de los puntos de interes  
          text(st_coordinates(p_int[[1]]), labels = p_int[[1]]$nombre, pos = 4,family="Old English Text MT",col = 'white',cex=0.9)
        text(st_coordinates(p_int[[2]]), labels = p_int[[2]]$nombre, pos = 4,family="Old English Text MT",col = 'white',cex=0.7)
      # añade escala y norte   
      addscalebar()#(style="ticks")
      addnortharrow(pos = "topright", scale = 0.5)
    }    

Resultado

Juntamos el proceso en otra función que llamaremos mapa del tesoro y que nos genera el mapa deseado, de manera aleatoria y en un santiamén.

Como es un blog personal, he decidido llamar a estos mapas los del tesoro del capitán Vilber, jeje:

################################################
# Funciones para hacer un mapa del tesoro con R
# Autor: Fernando Villalba Bergado
# Fecha: 2018-2019
################################################
    mapa_tesoro<-function(){
            radio<-10000
            #crea una isla principal
            islagrande<-crea_isla(R=radio)%>% pol_to_sf()%>%validar_contorno()
            # crea archipielago
            a<-archipielago(islagrande,radio=radio)
            #junta las islas
            islas1<-junta_igran_archi(islagrande,a)
            #da los nombres y los punto de interes
            nom_I<-nom_islas(length(a))
            pun_int<-p_interes(islagrande)
              
            #pinta_mapa(is2,nom_I)
            aux<-pinta_mapa_negro(islas1,nom_I,pun_int)
            return(aux)
          }

    
     mapa_tesoro()

## NULL

Cambiar los colores es bastante facil, aquí dejo una muestra con otras opcciones más coloridas:

## NULL

## NULL