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

21

Quiero trazar los límites de los países de América del Norte sobre una imagen de trama que representa alguna variable y luego superponer los contornos en la parte superior de la trama usando R. He tenido éxito al hacer esto utilizando gráficos de base y celosía, pero parece que proceso de trazado es demasiado lento! No he hecho esto en ggplot2 todavía, pero dudo que le vaya mejor en términos de velocidad.

Tengo los datos en un archivo netcdf creado a partir de un archivo de grib. Por ahora, descargué las fronteras del país para Canadá, Estados Unidos y México, que estaban disponibles en archivos RData de GADM que se lee en R como SpatialPolygonsDataFrame objetos.

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 forma de acelerar el trazado de los polígonos? En el sistema en el que estoy trabajando, el trazado lleva varios minutos. Quiero finalmente hacer una función que genere fácilmente un número de estas parcelas para la inspección, y creo que estaré trazando muchos de estos mapas, ¡así que quiero aumentar la velocidad de las parcelas!

¡Gracias!

    
pregunta ialm 30.05.2013 - 20:00

3 respuestas

27

Encontré 3 formas de aumentar la velocidad de trazado de los bordes del país desde archivos de formas para R. Encontré algo de inspiración y código de here y here .

(1) Podemos extraer las coordenadas de los archivos de formas para obtener la longitud y latitudes de los polígonos. Luego podemos colocarlos en un marco de datos con la primera columna que contiene las longitudes y la segunda columna que contiene latitudes. Las diferentes formas están separadas por NAs.

(2) Podemos eliminar algunos polígonos de nuestro archivo de formas. El archivo de formas es muy, muy detallado, pero algunas de las formas son pequeñas islas que no son importantes (para mis parcelas, de todos modos). Podemos establecer un umbral de área de polígono mínimo para mantener los polígonos más grandes.

(3) Podemos simplificar la geometría de nuestras formas mediante el uso de Douglas- Algoritmo de Peuker . Los bordes de nuestras formas poligonales pueden simplificarse, ya que son muy intrincados en el archivo original. Afortunadamente, hay un paquete, rgeos , que implementa esto.

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: extraiga las coordenadas de los archivos de formas en un marco de datos y dibuje líneas

La principal desventaja es que perdemos algo de información aquí en comparación con mantener el objeto como un objeto SpatialPolygonsDataFrame, como la proyección. Sin embargo, podemos volverlo a convertir 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 aproximadamente 2 millones de filas.

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([email protected])

  # 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([email protected][[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, 
                       [email protected][[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: eliminar polígonos pequeños

Hay muchas islas pequeñas que no son muy importantes. Si verifica 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 los polígonos para 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([email protected], 
                  function(x) sapply([email protected], function(y) [email protected]))

  # 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 and extract them
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]] >= 1){
      [email protected][[i]]@Polygons <- [email protected][[i]]@Polygons[bigpolys[[i]]]
      [email protected][[i]]@plotOrder <- 1:length([email protected][[i]]@Polygons)
    }
  }
  return(poly)
}

Método 3: simplifica la geometría de las formas de polígono

Podemos reducir el número de vértices en nuestras formas poligonales usando la función gSimplify del paquete rgeos

Código:

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

Algunos puntos de referencia:

Usé system.time transcurrido para comparar mis tiempos de trazado. Tenga en cuenta que estos son solo los tiempos para trazar los países, sin las curvas de nivel y otras cosas adicionales. Para los objetos sp, acabo de usar la función plot . Para los objetos del marco de datos, usé la función plot con type='l' y la función lines .

Trazando los polígonos originales de Canadá, EE. UU., México:

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 proporciona aumentos de velocidad suficientes para trazar polígonos. El uso de los métodos 2 + 3 + 1 agrega el problema de perder las propiedades agradables de los objetos sp , y mi principal dificultad es la aplicación de proyecciones. Hacké algo juntos para proyectar un objeto de marco de datos, pero se ejecuta bastante lento. Creo que usar el método 2 + 3 me proporciona suficientes aceleraciones hasta que puedo solucionar los problemas usando el método 2 + 3 + 1.

    
respondido por el ialm 31.05.2013 - 23:12
2

Todo el mundo debería considerar la transferencia a un paquete sf (características 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 lectura en un shp y trazado a través de ggplot2.

Nota: necesitas reinstalar ggplot2 desde la compilación más reciente en github (ver más abajo)

library(rgdal)
library(sp)
library(sf)
library(plyr)
devtools::install_github("tidyverse/ggplot2")
library(ggplot2)

# Load the shape files
can<-getData('GADM', country="CAN", level=0)
td <- file.path(tempdir(), "rgdal_examples"); dir.create(td)
st_write(st_as_sf(can),file.path(td,'can.shp'))


ptm <- proc.time()
  can = readOGR(dsn=td, layer="can")
  [email protected]$id = rownames([email protected])
  can.points = fortify(can, region="id")
  can.df = join(can.points, [email protected], by="id")
  ggplot(can.df) +  geom_polygon(aes(long,lat,group=group,fill='NAME_ENGLISH'))
proc.time() - ptm

user  system elapsed 
683.344   0.980 684.51 

ptm <- proc.time()
  can2 = st_read(file.path(td,'can.shp'))  
  ggplot(can2)+geom_sf( aes(fill = 'NAME_ENGLISH' )) 
proc.time() - ptm

user  system elapsed 
11.340   0.096  11.433 
    
respondido por el mmann1123 21.12.2017 - 00:20
0

Los datos GADM tienen una resolución espacial muy alta de las líneas costeras. Si no necesita, puede utilizar un conjunto de datos más generalizado. Los enfoques de ialm son muy interesantes, pero una alternativa simple es utilizar los datos 'wrld_simpl' que vienen con 'maptools'

library(maptools)
data(wrld_simpl)
plot(wrld_simpl)
    
respondido por el Robert Hijmans 01.06.2013 - 19:22

Lea otras preguntas en las etiquetas