Quiero trazar las fronteras de los países de América del Norte sobre una imagen de trama que representa algunas variables y luego superponer los contornos en la parte superior de la trama con R. He tenido éxito al hacer esto usando gráficos básicos y celosía, pero parece que el proceso de trazado es Demasiado lento! Todavía no he hecho esto en ggplot2, pero dudo que le vaya mejor en términos de velocidad.
Tengo los datos en un archivo netcdf creado a partir de un archivo grib. Por ahora, descargué las fronteras del país para Canadá, EE. UU. Y México, que estaban disponibles en archivos RData de GADM que se leían en R como objetos SpatialPolygonsDataFrame.
Aquí hay un código:
# Load packages
library(raster)
#library(ncdf) # If you cannot install ncdf4
library(ncdf4)
# Read in the file, get the 13th layer
# fn <- 'path_to_file'
r <- raster(fn, band=13)
# Set the projection and extent
p4 <- "+proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +x_0=32.46341 +y_0=32.46341 +lon_0=-107 +lat_0=1.0"
projection(r) <- CRS(p4)
extent(r) <- c(-5648.71, 5680.72, 1481.40, 10430.62)
# Get the country borders
# This will download the RData files to your working directory
can<-getData('GADM', country="CAN", level=1)
usa<-getData('GADM', country="USA", level=1)
mex<-getData('GADM', country="MEX", level=1)
# Project to model grid
can_p <- spTransform(can, CRS(p4))
usa_p <- spTransform(usa, CRS(p4))
mex_p <- spTransform(mex, CRS(p4))
### USING BASE GRAPHICS
par(mar=c(0,0,0,0))
# Plot the raster
bins <- 100
plot(r, axes=FALSE, box=FALSE, legend=FALSE,
col=rev( rainbow(bins,start=0,end=1) ),
breaks=seq(4500,6000,length.out=bins))
plot(r, legend.only=TRUE, col=rev( rainbow(bins,start=0,end=1)),
legend.width=0.5, legend.shrink=0.75,
breaks=seq(4500,6000,length.out=bins),
axis.args=list(at=seq(4500,6000,length.out=11),
labels=seq(4500,6000,length.out=11),
cex.axis=0.5),
legend.args=list(text='Height (m)', side=4, font=2,
line=2, cex=0.8))
# Plot the borders
# These are so slow!!
plot(can_p, add=TRUE, border='white', lwd=2)
plot(usa_p, add=TRUE, border='white', lwd=2)
plot(mex_p, add=TRUE, border='white', lwd=2)
# Add the contours
contour(r, add=TRUE, nlevel=5)
### USING LATTICE
library(rasterVis)
# Some settings for our themes
myTheme <- RdBuTheme()
myTheme$axis.line$col<-"transparent"
myTheme$add.line$alpha <- 1
myTheme2 <- myTheme
myTheme2$regions$col <- 'transparent'
myTheme2$add.text$cex <- 0.7
myTheme2$add.line$lwd <- 1
myTheme2$add.line$alpha <- 0.8
# Get JUST the contour lines
contours <- contourplot(r, margin=FALSE, scales=list(draw=FALSE),
par.settings=myTheme2, pretty=TRUE, key=NULL, cuts=5,
labels=TRUE)
# Plot the colour
levels <- levelplot(r, contour=FALSE, margin=FALSE, scales=list(draw=FALSE),
par.settings = myTheme, cuts=100)
# Plot!
levels +
layer(sp.polygons(can_p, col='green', lwd=2)) +
layer(sp.polygons(usa_p, col='green', lwd=2)) +
layer(sp.polygons(mex_p, col='green', lwd=2)) +
contours
¿Hay alguna manera de acelerar el trazado de los polígonos? En el sistema en el que estoy trabajando, el trazado lleva varios minutos. Finalmente, quiero hacer una función que genere fácilmente una serie de estas parcelas para inspección, y creo que trazaré muchos de estos mapas, ¡así que quiero aumentar la velocidad de las parcelas!
¡Gracias!
Respuestas:
Encontré 3 formas de aumentar la velocidad de trazar las fronteras del país a partir de archivos de formas para R. Encontré algo de inspiración y código aquí y aquí .
(1) Podemos extraer las coordenadas de los archivos de formas para obtener la longitud y latitudes de los polígonos. Luego podemos ponerlos en un marco de datos con la primera columna que contiene las longitudes y la segunda columna que contiene las latitudes. Las diferentes formas están separadas por NA.
(2) Podemos eliminar algunos polígonos de nuestro archivo de forma. El archivo de formas es muy, muy detallado, pero algunas de las formas son pequeñas islas que no son importantes (para mis tramas, de todos modos). Podemos establecer un umbral mínimo de área de polígono para mantener los polígonos más grandes.
(3) Podemos simplificar la geometría de nuestras formas usando el algoritmo Douglas-Peuker . Los bordes de nuestras formas poligonales se pueden simplificar, ya que son muy intrincados en el archivo original. Afortunadamente, hay un paquete
rgeos
que implementa esto.Preparar:
Método 1: Extraiga las coordenadas de los archivos de forma en un marco de datos y líneas de trazado
La principal desventaja es que perdemos cierta información aquí en comparación con mantener el objeto como un objeto SpatialPolygonsDataFrame, como la proyección. Sin embargo, podemos volver a convertirlo en un objeto sp y volver a agregar la información de proyección, y aún es más rápido que trazar los datos originales.
Tenga en cuenta que este código se ejecuta muy lentamente en el archivo original porque hay muchas formas y el marco de datos resultante tiene una longitud de ~ 2 millones de filas.
Código:
Método 2: eliminar los polígonos pequeños
Hay muchas islas pequeñas que no son muy importantes. Si revisa algunos de los cuantiles de las áreas para los polígonos, vemos que muchos de ellos son minúsculos. Para la trama de Canadá, pasé de trazar más de mil polígonos a solo cientos de polígonos.
Cuantiles para el tamaño de polígonos para Canadá:
Código:
Método 3: simplifique la geometría de las formas poligonales
Podemos reducir el número de vértices en nuestras formas poligonales usando la
gSimplify
función delrgeos
paqueteCódigo:
Algunos puntos de referencia:
Solía transcurrir
system.time
para comparar mis tiempos de trazado. Tenga en cuenta que estos son solo los momentos para trazar los países, sin las líneas de contorno y otras cosas adicionales. Para los objetos sp, acabo de usar laplot
función. Para los objetos del marco de datos, utilicé laplot
función contype='l'
y lalines
función.Trazando los polígonos originales de Canadá, Estados Unidos y México:
73,009 segundos
Usando el Método 1:
2.449 segundos
Usando el Método 2:
17.660 segundos
Usando el Método 3:
16.695 segundos
Usando el Método 2 + 1:
1.729 segundos
Usando el Método 2 + 3:
0.445 segundos
Usando el Método 2 + 3 + 1:
0,172 segundos
Otras observaciones:
Parece que la combinación de los métodos 2 + 3 proporciona aceleraciones suficientes para el trazado de polígonos. Usar los métodos 2 + 3 + 1 agrega el problema de perder las bonitas propiedades de los
sp
objetos, y mi principal dificultad es aplicar proyecciones. Pirateé algo juntos para proyectar un objeto de marco de datos, pero funciona bastante lento. Creo que usar el método 2 + 3 proporciona suficientes aceleraciones para mí hasta que pueda obtener los inconvenientes de usar el método 2 + 3 + 1.fuente
Todos deberían considerar la transferencia al paquete sf (funciones espaciales) en lugar de sp. Es significativamente más rápido (1/60 en este caso) y más fácil de usar. Aquí hay un ejemplo de leer en un shp y trazar a través de ggplot2.
Nota: debe reinstalar ggplot2 desde la compilación más reciente en github (consulte a continuación)
fuente
Los datos de GADM tienen una resolución espacial muy alta de las costas. Si no lo necesita, puede usar un conjunto de datos más generalizado. Los enfoques de ialm son muy interesantes, pero una alternativa simple es usar los datos 'wrld_simpl' que vienen con 'maptools'
fuente