Cómo acelerar el trazado de polígonos en R?

Question

Quiero trazar las fronteras de los países de América del Norte a través de una imagen ráster que representa una variable y, a continuación, superposición de contornos en la parte superior de la parcela, mediante R. he tenido éxito en hacer esto usando la base de los gráficos y la red, pero parece que el trazado proceso es demasiado lento! No he hecho esto en ggplot2, pero dudo que se les va mejor en términos de velocidad.

Tengo los datos en un archivo netcdf creado a partir de un archivo grib. Por ahora, he descargado las fronteras de los países de Canadá, Estados Unidos y México, que estaban disponibles en RData archivos de GADM que leer en R como SpatialPolygonsDataFrame objetos.

Aquí está el 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

Es allí una manera de acelerar el trazado de los polígonos? En el sistema que estoy trabajando, que la conspiración lleva varios minutos. Quiero, finalmente, hacer una función que va a generar fácilmente un número de estas parcelas para la inspección, y creo que me va a estar conspirando muchos de estos mapas, por lo que quiero aumentar la velocidad de las parcelas!

Gracias!

Answer 1

He encontrado 3 formas de aumentar la velocidad de trazado de las fronteras del país de archivos shape de R. he encontrado un poco de inspiración y de código a partir de aquí y aquí.

(1) se pueden extraer las coordenadas de la forma de archivos para obtener la longitud y latitudes de los polígonos. A continuación, los podemos poner en un marco de datos en la primera columna contiene las longitudes y la segunda columna contiene las latitudes. Las diferentes formas están separados por NAs.

(2) podemos eliminar algunos polígonos de nuestro archivo de forma. El archivo de forma es muy, muy detallado, pero algunas de las formas son pequeñas islas que no son importantes (para mi maquina, de todos modos). Podemos establecer un mínimo polígono de la zona de umbral para mantener el más grande de los polígonos.

(3) se puede simplificar la geometría de nuestras formas mediante las Douglas-Peuker algoritmo. Los bordes de nuestro polígono formas puede ser simplificado, ya que son muy complicados en el archivo original. Afortunadamente, no es un paquete, rgeos, que implementa este.

Configurar:

# Load packages
library(rgdal)
library(raster)
library(sp)
library(rgeos)

# Load the shape files
can<-getData('GADM', country="CAN", level=0)
usa<-getData('GADM', country="USA", level=0)
mex<-getData('GADM', country="MEX", level=0)

Método 1: Extraer las coordenadas de la forma de archivos en un marco de datos y líneas argumentales

La principal desventaja es que perdemos aquí algo de información cuando se compara a mantener el objeto como un SpatialPolygonsDataFrame objeto, tales como la proyección. Sin embargo, podemos convertirlo de nuevo en un sp objeto y volver a agregar la información de proyección, y es aún más rápido que el trazado de la original de datos.

Tenga en cuenta que este código se ejecuta muy lentamente en el archivo original, porque hay un montón de formas, y el resultado de la trama de datos es de ~2 millones de líneas.

Código:

# Convert the polygons into data frames so we can make lines
poly2df <- function(poly) {
  # Convert the polygons into data frames so we can make lines
  # Number of regions
  n_regions <- length(poly@polygons)

  # Get the coords into a data frame
  poly_df <- c()
  for(i in 1:n_regions) {
    # Number of polygons for first region
    n_poly <- length(poly@polygons[[i]]@Polygons)
    print(paste("There are",n_poly,"polygons"))
    # Create progress bar
    pb <- txtProgressBar(min = 0, max = n_poly, style = 3)
    for(j in 1:n_poly) {
      poly_df <- rbind(poly_df, NA, 
                       poly@polygons[[i]]@Polygons[[j]]@coords)
      # Update progress bar
      setTxtProgressBar(pb, j)
    }
    close(pb)
    print(paste("Finished region",i,"of",n_regions))
  }
  poly_df <- data.frame(poly_df)
  names(poly_df) <- c('lon','lat')
  return(poly_df)
}

Método 2: Quitar los pequeños polígonos

Hay muchas pequeñas islas que no son muy importantes. Si usted comprobar algunos de los cuantiles de las áreas de los polígonos, vemos que muchos de ellos son minúsculos. Para el Canadá de la trama, me fui abajo de conspirar más de un millar de polígonos a cientos de polígonos.

Cuantiles para el tamaño de los polígonos de Canadá:

          0%          25%          50%          75%         100% 
4.335000e-10 8.780845e-06 2.666822e-05 1.800103e-04 2.104909e+02 

Código:

# Get the main polygons, will determine by area.
getSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons, 
                  function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][4] >= 1){
      poly@polygons[[i]]@Polygons <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:length(poly@polygons[[i]]@Polygons)
    }
  }
  return(poly)
}

Método 3: Simplificar la geometría del polígono formas

Podemos reducir el número de vértices en nuestro polígono formas mediante las gSimplify función de la rgeos paquete

Código:

can <- getData('GADM', country="CAN", level=0)
can <- gSimplify(can, tol=0.01, topologyPreserve=TRUE)

Algunos puntos de referencia:

He utilizado transcurrido system.time de punto de referencia para mi conspirar veces. Tenga en cuenta que estos son sólo los tiempos para el trazado de los países, sin las líneas de contorno y otras cosas. Para el sp objetos, sólo utiliza la plot función. Para el marco de datos de objetos, he utilizado el plot función de con type='l' y lines función.

El trazado de la original de Canadá, Estados Unidos, México polígonos:

73.009 segundos

Utilizando El Método 1:

2.449 segundos

Utilizando El Método 2:

17.660 segundos

Utilizando El Método 3:

16.695 segundos

Utilizando El Método 2 + 1:

1.729 segundos

Utilizando El Método 2 + 3:

0.445 segundos

Utilizando El Método 2 + 3 + 1:

0.172 segundos

Otras observaciones:

Parece que la combinación de los métodos 2 + 3 da la suficiente velocidad de ups para el trazado de los polígonos. El uso de métodos 2 + 3 + 1 añade el problema de la pérdida de las propiedades atractivas de sp objetos, y mi principal dificultad es la aplicación de las proyecciones. He hackeado algo juntos para el proyecto de un marco de datos de objeto, pero corre bastante lento. Creo que utilizando el método 2 + 3 proporciona la suficiente velocidad de ups para mí hasta que yo pueda obtener los kinks a cabo utilizando el método de 2 + 3 + 1.

Answer 2

GADM de datos tienen una muy alta resolución espacial de las costas. Si usted no necesita que usted puede utilizar una forma más generalizada del conjunto de datos. ialm enfoques son muy interesantes, pero una alternativa sencilla es usar el wrld_simpl de datos que viene con 'maptools'

library(maptools)
data(wrld_simpl)
plot(wrld_simpl)