Bagaimana mempercepat plot poligon di R?

Saya ingin memplot batas-batas negara Amerika Utara pada gambar raster yang menggambarkan beberapa variabel dan kemudian overlay kontur di atas plot menggunakan R. Saya telah berhasil melakukan ini dengan menggunakan grafik dasar dan kisi, tetapi tampaknya proses ploting adalah terlalu lambat! Saya belum melakukan ini di ggplot2, tapi saya ragu itu akan lebih baik dalam hal kecepatan.

Saya memiliki data dalam file netcdf yang dibuat dari file grib. Untuk saat ini, saya mengunduh perbatasan negara untuk Kanada, Amerika Serikat, dan Meksiko, yang tersedia dalam file RData dari GADM yang berbunyi R sebagai objek SpatialPolygonsDataFrame.

Ini beberapa kode:

# 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

Apakah ada cara untuk mempercepat perencanaan poligon? Pada sistem yang sedang saya kerjakan, perencanaan membutuhkan waktu beberapa menit. Saya ingin akhirnya membuat fungsi yang akan dengan mudah menghasilkan sejumlah plot ini untuk diperiksa, dan saya rasa saya akan memplot banyak peta ini, jadi saya ingin meningkatkan kecepatan plot!

Terima kasih!

Saya menemukan 3 cara untuk meningkatkan kecepatan merencanakan batas negara dari file bentuk untuk R. Saya menemukan beberapa inspirasi dan kode dari sini dan di sini .

(1) Kita dapat mengekstrak koordinat dari file bentuk untuk mendapatkan garis bujur dan garis lintang poligon. Kemudian kita bisa meletakkannya dalam bingkai data dengan kolom pertama berisi garis bujur dan kolom kedua berisi garis lintang. Bentuk yang berbeda dipisahkan oleh NAS.

(2) Kami dapat menghapus beberapa poligon dari file bentuk kami. File bentuk sangat, sangat rinci, tetapi beberapa bentuknya adalah pulau-pulau kecil yang tidak penting (untuk plot saya, lagipula). Kita dapat menetapkan ambang area poligon minimum untuk menjaga poligon yang lebih besar.

(3) Kita dapat menyederhanakan geometri bentuk kita menggunakan algoritma Douglas-Peuker . Tepi bentuk poligon kami dapat disederhanakan, karena sangat rumit dalam file aslinya. Untungnya, ada paket rgeos,, yang mengimplementasikan ini.

Mendirikan:

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

Metode 1: Ekstrak koordinat dari file bentuk ke dalam bingkai data dan garis petak

Kerugian utama adalah bahwa kami kehilangan beberapa informasi di sini jika dibandingkan dengan menjaga objek sebagai objek SpatialPolygonsDataFrame, seperti proyeksi. Namun, kita dapat mengubahnya kembali menjadi objek sp dan menambahkan kembali informasi proyeksi, dan itu masih lebih cepat daripada memplot data asli.

Perhatikan bahwa kode ini berjalan sangat lambat pada file asli karena ada banyak bentuk, dan bingkai data yang dihasilkan ~ 2 juta baris.

Kode:

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

Metode 2: Hapus poligon kecil

Ada banyak pulau kecil yang tidak terlalu penting. Jika Anda memeriksa beberapa kuantil area untuk poligon, kami melihat bahwa banyak di antaranya sangat kecil. Untuk plot Kanada, saya turun dari memplot lebih dari seribu poligon menjadi hanya ratusan poligon.

Jumlah untuk ukuran poligon untuk Kanada:

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

Kode:

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

Metode 3: Sederhanakan geometri bentuk poligon

Kita dapat mengurangi jumlah simpul dalam bentuk poligon kami menggunakan gSimplifyfungsi dari rgeospaket

Kode:

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

Beberapa tolok ukur:

Saya menggunakan berlalu system.timeuntuk tolok ukur waktu merencanakan saya. Perhatikan bahwa ini hanya waktu untuk merencanakan negara, tanpa garis kontur dan hal-hal tambahan lainnya. Untuk objek sp, saya hanya menggunakan plotfungsinya. Untuk objek bingkai data, saya menggunakan plotfungsi type='l'dan linesfungsinya.

Merencanakan poligon Kanada, Amerika Serikat, Meksiko asli:

73,009 detik

Menggunakan Metode 1:

2,449 detik

Menggunakan Metode 2:

17,660 detik

Menggunakan Metode 3:

16.695 detik

Menggunakan Metode 2 + 1:

1,729 detik

Menggunakan Metode 2 + 3:

0,445 detik

Menggunakan Metode 2 + 3 + 1:

0,172 detik

Keterangan lain:

Tampaknya kombinasi metode 2 + 3 memberikan kecepatan yang cukup untuk memplot poligon. Menggunakan metode 2 + 3 + 1 menambah masalah kehilangan properti bagus spobjek, dan kesulitan utama saya adalah menerapkan proyeksi. Saya meretas sesuatu bersama-sama untuk memproyeksikan objek bingkai data, tetapi itu berjalan agak lambat. Saya pikir menggunakan metode 2 + 3 memberikan kecepatan yang cukup bagi saya sampai saya bisa mendapatkan ketegaran menggunakan metode 2 + 3 + 1.

Semua orang harus melihat ke transfer ke paket sf (fitur spasial) alih-alih sp. Secara signifikan lebih cepat (1/60 dalam hal ini) dan lebih mudah digunakan. Berikut adalah contoh membaca dalam shp dan memplot via ggplot2.

Catatan: Anda perlu menginstal ulang ggplot2 dari build terbaru di github (lihat di bawah)

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 

Data GADM memiliki resolusi spasial yang sangat tinggi dari garis pantai. Jika tidak perlu, Anda dapat menggunakan kumpulan data yang lebih umum. Pendekatan ialm sangat menarik, tetapi alternatif sederhana adalah dengan menggunakan data 'wrld_simpl' yang dilengkapi dengan 'maptools'

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