A continuación se muestra el código R utilizado para crear mapas 3D con datos obtenidos de la Paleobiology Database para explorar patrones de extinción
Primero, descargaremos los datos que vamos a graficar. Para este ejemplo usaremos los géneros del orden Proboscidea. Los datos son filtrados para remover las columnas que no nos son útiles. Todos aquellos registros sin un género asignado son removidos también.
#Download dataset from Paleobiology Database. The id in the url corresponds to the taxa ID in the PBDB
url <- 'https://paleobiodb.org/data1.2/occs/list.txt?base_id=43230&show=full'
full_data <- read_csv(file = url)
sumarized <- full_data [,c(25,14,16,17,30,31)]
#Filter data excluding occurrences without genus
complete <- subset(sumarized, (!is.na(sumarized$genus))&!is.na(sumarized$lat))
Es necesario crear un mapa base para usarse como fondo. En este caso, se usará un mapa global, por lo que no es necesario ajustar los límites del mapa. Es necesario remover todos los elementos como las rejillas para generar un mapa más limpio.
A continuación es necesario definir los límites X y Y en el mapa. Para posteriormente exportar el mapa como una imagen PNG.
#Define limits of the map
xlim = ggplot_build(basemap)$layout$panel_scales_x[[1]]$range$range
ylim = ggplot_build(basemap)$layout$panel_scales_y[[1]]$range$range
#Save map as a PNG object
ggsave("map_img.png",plot = basemap, width = diff(xlim), height = diff(ylim), units = "mm")
#Load the PNG map object
background <- readPNG("map_img.png")
Basemap
El mapa ahora puede ser graficado con todos los puntos de ocurrencias sobre él. Diferentes formas y tamaños deben ser asignados a los valores de mínima y máxima edad. El valor de edad mínima debe ser más pequeño para una mejor visualización
#Map with points. Different shape and size for max_ma and min_ma are necesary for a better visualization
map1 <- ggplot(data = region) + annotation_custom(rasterGrob(background, width=unit(1,"npc"), height=unit(1,"npc")), -Inf, Inf, -Inf, Inf)+xlim(xlim[1],xlim[2]) + ylim(ylim[1],ylim[2])+
geom_point(data = complete, aes(x = lng, y = lat, colour=max_ma),shape=5, size=0.2) +
geom_point(data = complete, aes(x = lng, y = lat, colour=min_ma),shape=3, size=1)+
scale_colour_gradient(name = "Millions of years",
limits=range(complete$max_ma))+
ggtitle("occurreces")+
theme(plot.title = element_text(hjust = 0.5))
Para guardar el mapa 3D usaremos el siguiente comando SIN CERRAR LA VENTANA donde nos muestra la vista previa.
#Save 3D object
save_obj("plot3D.obj", save_texture = TRUE)
Usando más código y opciones de personalización, es posible elaborar mapas más complejos. En el siguiente ejemplo, se muestra un mapa de proboscideos por país.
url <- 'https://paleobiodb.org/data1.2/occs/list.txt?base_id=105150&show=full'
full_data <- read_csv(file = url)
sumarized <- full_data [,c(6,7,10,12,14,15,16,17,24,25,30,31,33,34,35,36,37,38,43,44,45,51)]
complete <- subset(sumarized, (!is.na(sumarized$family))&(!sumarized$family=="NO_FAMILY_SPECIFIED"))
names(complete)[16] <- "id"
#Due to political changes, all occureces of Western Sahara, will be considered as Morocco occureces
complete$id[complete$id =="EH"] <- "MA"
#Ocurece by country
fam_by_country <- aggregate(data = complete,
family ~ id,
function(x) length(unique(x)))
occ_by_country <- split(complete, complete[["id"]])
#Unique values
genus_unique <-unique(complete[c("genus")])
genus_unique <- na.omit(genus_unique)
family_unique <-unique(complete[c("family")])
family_unique <-na.omit(family_unique)
#Map
world <- ne_countries(scale = "medium", returnclass = "sf")
names(world)[44] <- "id"
fam_by_country$id[fam_by_country$id =="UK"] <- "GB"
world <- merge(world,fam_by_country, by="id", all=TRUE)
world$family[is.na(world$family)]<- 0
#Match dataset with centroid country
wmap <- getMap(resolution="high")
centroids <- gCentroid(wmap, byid=TRUE)
df <- as.data.frame(centroids)
df <- tibble::rownames_to_column(df, "country")
df$id <-countrycode::countrycode(sourcevar = df$country, origin = "country.name", destination ="iso2c", custom_match = c("United Kingdom"="UK"))
df <- subset(df, !is.na(df$id)&(!df$country== "US Naval Base Guantanamo Bay" & !df$country== "Korea No Mans Area" & !df$country== "Cyprus No Mans Area" & !df$country=="Northern Cyprus" & !df$country=="West Bank"))
merged <- merge(df,fam_by_country, by="id")
merged$id[merged$id =="UK"] <- "GB"
#convert map to png
map <- ggplot(data = world) +
geom_sf(aes(fill = family), colour="#FFFFFF") +
scale_fill_gradient(limits=range(world$family),low="#FFF3B0", high="#E09F3E") +
theme(legend.position = "none",axis.line=element_blank(),axis.text.x=element_blank(), axis.title.x=element_blank(),axis.text.y=element_blank(), axis.title.y=element_blank(),axis.ticks=element_blank(),panel.background = element_blank())
xlim = ggplot_build(map)$layout$panel_scales_x[[1]]$range$range
ylim = ggplot_build(map)$layout$panel_scales_y[[1]]$range$range
ggsave("map_img.png",plot = map, width = diff(xlim), height = diff(ylim), units = "mm")
worldmap <- readPNG("map_img.png")
mapa1 <- ggplot(merged) + annotation_custom(rasterGrob(worldmap, width=unit(1,"npc"), height=unit(1,"npc")), -Inf, Inf, -Inf, Inf) +
xlim(xlim[1],xlim[2]) + ylim(ylim[1],ylim[2]) +
geom_point(data = merged, aes(x = x, y = y, colour=family), size=2)+
scale_colour_gradient(name = "Number of families",
limits=range(merged$family),low="#FCB9B2", high="#B23A48") +
ggtitle("Mammuthidae family members per country")+
theme(legend.position = "bottom",axis.line=element_blank(),
axis.text.x=element_blank(), axis.title.x=element_blank(),
axis.text.y=element_blank(), axis.title.y=element_blank(),
axis.ticks=element_blank(),panel.background = element_blank(),plot.title = element_text(hjust = 0.5))
#3D Map
map_3d <- plot_gg(mapa1,multicore = TRUE, width = diff(xlim)/30 ,height=diff(ylim)/30, raytrace = FALSE, windowsize = c(1200, 960), fov = 70, zoom = 0.25, theta = 330, phi = 40)
Miguel Díaz de León
Lic. en Biología por la Universidad de Guadalajara (2010-2015)
Maestría en Ciencias de La Tierra UNAM (2017-2020)
