Wie beschleunige ich das Zeichnen von Polygonen in R?

Ich möchte die Landesgrenzen Nordamerikas über ein Rasterbild zeichnen, das eine Variable darstellt, und dann Konturen mit R über das Diagramm legen. Es ist mir gelungen, dies unter Verwendung von Basisgrafiken und Gittern zu tun, aber es scheint, dass dies der Vorgang ist viel zu langsam! Ich habe dies in ggplot2 noch nicht getan, bezweifle aber, dass es in Bezug auf die Geschwindigkeit besser abschneiden wird.

Ich habe die Daten in einer Netcdf-Datei aus einer Grib-Datei erstellt. Im Moment habe ich die Ländergrenzen für Kanada, USA und Mexiko heruntergeladen, die in RData-Dateien von GADM verfügbar waren, die als SpatialPolygonsDataFrame-Objekte in R eingelesen wurden.

Hier ist ein Code:

# 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

Gibt es eine Möglichkeit, das Zeichnen der Polygone zu beschleunigen? Auf dem System, an dem ich arbeite, dauert das Plotten einige Minuten. Ich möchte irgendwann eine Funktion erstellen, mit der eine Reihe dieser Diagramme einfach zur Überprüfung erstellt werden kann, und ich gehe davon aus, dass ich viele dieser Karten zeichnen werde, also möchte ich die Geschwindigkeit der Diagramme erhöhen!

Vielen Dank!

Ich habe drei Möglichkeiten gefunden, um die Geschwindigkeit des Zeichnens der Ländergrenzen aus Formdateien für R zu erhöhen. Ich habe hier und hier Inspiration und Code gefunden .

(1) Wir können die Koordinaten aus den Formdateien extrahieren, um die Längen- und Breitengrade der Polygone zu erhalten. Dann können wir sie in einen Datenrahmen einfügen, wobei die erste Spalte die Längengrade und die zweite Spalte die Breitengrade enthält. Die verschiedenen Formen werden durch NAs getrennt.

(2) Wir können einige Polygone aus unserer Formdatei entfernen. Die Formdatei ist sehr, sehr detailliert, aber einige der Formen sind winzige Inseln, die unwichtig sind (für meine Zeichnungen jedenfalls). Wir können einen Mindestschwellenwert für die Polygonfläche festlegen, um die größeren Polygone beizubehalten.

(3) Mit dem Douglas-Peuker-Algorithmus können wir die Geometrie unserer Formen vereinfachen . Die Kanten unserer Polygonformen können vereinfacht werden, da sie in der Originaldatei sehr kompliziert sind. Zum Glück gibt es ein Paket rgeos, das dies implementiert.

Installieren:

# 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)

Methode 1: Extrahieren Sie die Koordinaten aus den Formdateien in einen Datenrahmen und zeichnen Sie Linien

Der Hauptnachteil besteht darin, dass wir hier einige Informationen verlieren, wenn wir das Objekt als SpatialPolygonsDataFrame-Objekt behalten, z. B. die Projektion. Wir können es jedoch wieder in ein sp-Objekt umwandeln und die Projektionsinformationen hinzufügen, und es ist immer noch schneller als das Plotten der Originaldaten.

Beachten Sie, dass dieser Code in der Originaldatei sehr langsam ausgeführt wird, da viele Formen vorhanden sind und der resultierende Datenrahmen ~ 2 Millionen Zeilen lang ist.

Code:

# 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)
}

Methode 2: Entfernen Sie kleine Polygone

Es gibt viele kleine Inseln, die nicht sehr wichtig sind. Wenn Sie einige der Quantile der Bereiche auf Polygone überprüfen, sehen wir, dass viele von ihnen winzig sind. Für den Kanada-Plot ging ich vom Plotten von über tausend Polygonen auf nur Hunderte von Polygonen zurück.

Quantile für die Größe von Polygonen für Kanada:

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

Code:

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

Methode 3: vereinfachen Sie die Geometrie der Polygonformen

Mit der gSimplifyFunktion aus dem rgeosPaket können wir die Anzahl der Eckpunkte in unseren Polygonformen reduzieren

Code:

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

Einige Benchmarks:

Ich habe die abgelaufenen system.timeZeiten verwendet, um meine Planungszeiten zu bestimmen. Beachten Sie, dass dies nur die Zeiten für die Darstellung der Länder sind, ohne die Konturlinien und andere zusätzliche Dinge. Für die sp-Objekte habe ich nur die plotFunktion verwendet. Für die Datenrahmenobjekte habe ich die plotFunktion mit type='l'und die linesFunktion verwendet.

Zeichnen der ursprünglichen Kanada-, USA- und Mexiko-Polygone:

73.009 Sekunden

Mit Methode 1:

2,449 Sekunden

Mit Methode 2:

17,660 Sekunden

Mit Methode 3:

16,695 Sekunden

Mit Methode 2 + 1:

1,729 Sekunden

Mit Methode 2 + 3:

0,445 Sekunden

Mit Methode 2 + 3 + 1:

0,172 Sekunden

Sonstige Anmerkungen:

Es scheint, dass die Kombination der Methoden 2 + 3 die Darstellung von Polygonen ausreichend beschleunigt. Die Verwendung der Methoden 2 + 3 + 1 fügt das Problem hinzu, die schönen Eigenschaften von spObjekten zu verlieren , und meine Hauptschwierigkeit besteht darin, Projektionen anzuwenden. Ich habe etwas zusammen gehackt, um ein Datenrahmenobjekt zu projizieren, aber es läuft ziemlich langsam. Ich denke, dass die Verwendung von Methode 2 + 3 eine ausreichende Beschleunigung für mich darstellt, bis ich die Nachteile der Verwendung von Methode 2 + 3 + 1 beseitigen kann.

Jeder sollte die Umstellung auf das Paket sf (Spatial Features) anstelle von sp in Betracht ziehen. Es ist deutlich schneller (in diesem Fall 1/60) und einfacher zu bedienen. Hier ist ein Beispiel für das Einlesen eines shp und das Plotten über ggplot2.

Hinweis: Sie müssen ggplot2 aus der neuesten Version von github neu installieren (siehe unten).

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")
  can@data$id = rownames(can@data)
  can.points = fortify(can, region="id")
  can.df = join(can.points, can@data, 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 

GADM-Daten haben eine sehr hohe räumliche Auflösung der Küsten. Wenn Sie dies nicht benötigen, können Sie einen allgemeineren Datensatz verwenden. Die Ansätze von ialm sind sehr interessant, aber eine einfache Alternative ist die Verwendung der 'wrld_simpl'-Daten, die mit' maptools 'geliefert werden.

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