Jak przyspieszyć wykreślanie wielokątów w R?

Chcę wykreślić granice kraju Ameryki Północnej na obrazie rastrowym przedstawiającym pewną zmienną, a następnie nałożyć wykres na wierzch wykresu za pomocą R. Udało mi się to zrobić przy użyciu podstawowej grafiki i sieci, ale wydaje się, że proces kreślenia jest o wiele za wolno! Nie zrobiłem tego jeszcze w ggplot2, ale wątpię, czy będzie lepiej, jeśli chodzi o szybkość.

Mam dane w pliku netcdf utworzonym z pliku grib. Na razie pobrałem granice kraju dla Kanady, USA i Meksyku, które były dostępne w plikach RData z GADM, które wczytują R jako obiekty SpatialPolygonsDataFrame.

Oto kod:

# 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

Czy istnieje sposób na przyspieszenie kreślenia wielokątów? W systemie, nad którym pracuję, kreślenie zajmuje kilka minut. Chcę w końcu stworzyć funkcję, która z łatwością wygeneruje pewną liczbę tych wykresów do kontroli, i sądzę, że wykreślę wiele z tych map, więc chcę zwiększyć prędkość wykresów!

Dzięki!

Znalazłem 3 sposoby na zwiększenie szybkości kreślenia granic kraju z plików kształtów dla R. Znalazłem inspirację i kod tutaj i tutaj .

(1) Możemy wyodrębnić współrzędne z plików kształtów, aby uzyskać długość i szerokość geograficzną wielokątów. Następnie możemy umieścić je w ramce danych, w której pierwsza kolumna zawiera długości geograficzne, a druga kolumna zawiera szerokości geograficzne. Różne kształty są oddzielone przez NA.

(2) Możemy usunąć niektóre wielokąty z naszego pliku kształtu. Plik kształtów jest bardzo, bardzo szczegółowy, ale niektóre kształty to małe wyspy, które nie są ważne (w każdym razie dla moich wykresów). Możemy ustawić minimalny próg powierzchni wielokąta, aby zachować większe wielokąty.

(3) Możemy uprościć geometrię naszych kształtów za pomocą algorytmu Douglas-Peuker . Krawędzie naszych kształtów wielokątów można uprościć, ponieważ są one bardzo skomplikowane w oryginalnym pliku. Na szczęście istnieje pakiet rgeos, który to implementuje.

Ustawiać:

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

Metoda 1: Wyodrębnij współrzędne z plików kształtów do ramki danych i narysuj linie

Główną wadą jest to, że tracimy tutaj trochę informacji w porównaniu z utrzymywaniem obiektu jako obiektu SpatialPolygonsDataFrame, takiego jak rzut. Możemy jednak przekształcić go z powrotem w obiekt sp i dodać informacje o rzutowaniu, a to wciąż jest szybsze niż drukowanie oryginalnych danych.

Zauważ, że ten kod działa bardzo wolno na oryginalnym pliku, ponieważ istnieje wiele kształtów, a wynikowa ramka danych ma długość ~ 2 milionów wierszy.

Kod:

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

Metoda 2: Usuń małe wielokąty

Istnieje wiele małych wysp, które nie są bardzo ważne. Jeśli sprawdzisz niektóre kwantyle obszarów dla wielokątów, zobaczymy, że wiele z nich jest maleńkich. W przypadku spisku w Kanadzie przeszedłem z wykreślania ponad tysiąca wielokątów do zaledwie setek wielokątów.

Kwantyle dla wielkości wielokątów dla Kanady:

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

Kod:

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

Metoda 3: Uprość geometrię kształtów wielokąta

Możemy zmniejszyć liczbę wierzchołków w naszych kształtach wielokątów za pomocą gSimplifyfunkcji z rgeospakietu

Kod:

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

Niektóre punkty odniesienia:

Upłynęło, system.timeby porównać moje czasy kreślenia. Pamiętaj, że to tylko czas na wykreślanie krajów, bez konturów i innych dodatkowych rzeczy. Dla obiektów sp użyłem właśnie tej plotfunkcji. W przypadku ramki danych obiektów, że stosuje się plotfunkcję z type='l'i linesfunkcji.

Rysowanie oryginalnych wielokątów w Kanadzie, USA i Meksyku:

73,009 sekund

Przy użyciu metody 1:

2.449 sekund

Przy użyciu metody 2:

17.660 sekund

Przy użyciu metody 3:

16,695 sekund

Przy użyciu metody 2 + 1:

1,729 sekundy

Przy użyciu metody 2 + 3:

0,445 sekundy

Przy użyciu metody 2 + 3 + 1:

0,172 sekundy

Inne uwagi:

Wydaje się, że połączenie metod 2 + 3 zapewnia wystarczające przyspieszenie kreślenia wielokątów. Korzystanie z metod 2 + 3 + 1 dodaje problem utraty ładnych właściwości spobiektów, a moją główną trudnością jest stosowanie projekcji. Zhakowałem coś razem, aby wyświetlić obiekt ramki danych, ale działa on dość wolno. Myślę, że użycie metody 2 + 3 zapewnia mi wystarczające przyspieszenie, dopóki nie będę w stanie załatwić problemu przy użyciu metody 2 + 3 + 1.

Każdy powinien rozważyć przeniesienie do pakietu sf (funkcje przestrzenne) zamiast sp. Jest znacznie szybszy (w tym przypadku 1/60) i łatwiejszy w użyciu. Oto przykład czytania w shp i kreślenia za pomocą ggplot2.

Uwaga: Musisz ponownie zainstalować ggplot2 z najnowszej wersji github (patrz poniżej)

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 

Dane GADM mają bardzo wysoką rozdzielczość przestrzenną linii brzegowych. Jeśli nie potrzebujesz, możesz użyć bardziej ogólnego zestawu danych. Podejścia ialm są bardzo interesujące, ale prostą alternatywą jest użycie danych „wrld_simpl” dostarczanych z „maptools”

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