Как ускорить прорисовку полигонов в R?

Я хочу нанести границы страны Северной Америки на растровое изображение, изображающее некоторую переменную, а затем наложить контуры на верхнюю часть графика, используя R. Я успешно справился с этим, используя базовую графику и решетку, но, похоже, процесс построения графика слишком медленно! Я еще не сделал этого в ggplot2, но я сомневаюсь, что это будет лучше с точки зрения скорости.

У меня есть данные в файле netcdf, созданном из файла grib. На данный момент я загрузил границы стран для Канады, США и Мексики, которые были доступны в файлах RData из GADM, которые считываются в R как объекты SpatialPolygonsDataFrame.

Вот некоторый код:

# 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

Есть ли способ ускорить прорисовку полигонов? В системе, над которой я работаю, построение графика занимает несколько минут. Я хочу в конечном итоге создать функцию, которая будет легко генерировать ряд этих графиков для проверки, и я рассчитываю, что буду строить многие из этих карт, поэтому я хочу увеличить скорость графиков!

Благодарность!

Я нашел 3 способа увеличить скорость построения границ страны из файлов формы для R. Я нашел вдохновение и код здесь и здесь .

(1) Мы можем извлечь координаты из файлов формы, чтобы получить долготу и широту полигонов. Затем мы можем поместить их в кадр данных с первым столбцом, содержащим долготы, и вторым столбцом, содержащим широты. Различные формы разделены NA.

(2) Мы можем удалить несколько полигонов из нашего файла формы. Файл формы очень, очень подробный, но некоторые формы - это крошечные островки, которые не важны (в любом случае, для моих графиков). Мы можем установить минимальный порог площади полигона, чтобы сохранить большие полигоны.

(3) Мы можем упростить геометрию наших фигур, используя алгоритм Дугласа-Пейкера . Края наших многоугольников можно упростить, поскольку они очень сложны в исходном файле. К счастью, есть пакет rgeos, который реализует это.

Настроить:

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

Метод 1: Извлечение координат из файлов формы в фрейм данных и построение линий

Основным недостатком является то, что мы теряем некоторую информацию здесь по сравнению с сохранением объекта в виде объекта SpatialPolygonsDataFrame, такого как проекция. Однако мы можем превратить его обратно в объект sp и добавить обратно информацию о проекции, и это все же быстрее, чем вывод исходных данных.

Обратите внимание, что этот код выполняется очень медленно в исходном файле, поскольку в нем много фигур, а длина результирующего фрейма данных составляет ~ 2 миллиона строк.

Код:

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

Способ 2: удалить небольшие полигоны

Есть много небольших островков, которые не очень важны. Если вы проверите некоторые из квантилей областей для полигонов, мы увидим, что многие из них являются крошечными. Что касается канадского сюжета, я перешел от построения более тысячи полигонов к сотням полигонов.

Квантили по размеру полигонов для Канады:

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

Код:

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

Способ 3: упрощение геометрии многоугольников

Мы можем уменьшить количество вершин в наших многоугольниках, используя gSimplifyфункцию из rgeosпакета

Код:

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

Некоторые тесты:

Я использовал прошедшее время system.timeдля сравнения моих графиков. Обратите внимание, что это как раз время для составления графика стран без контурных линий и других дополнительных вещей. Для объектов sp я просто использовал plotфункцию. Для объектов фрейма данных я использовал plotфункцию type='l'и linesфункцию.

Построение оригинальных полигонов Канады, США, Мексики:

73,009 секунды

Используя метод 1:

2,494 секунды

Используя метод 2:

17,660 секунд

Используя метод 3:

16,695 секунд

Используя метод 2 + 1:

1,729 секунды

Используя метод 2 + 3:

0,445 секунды

Используя метод 2 + 3 + 1:

0,172 секунды

Другие замечания:

Кажется, что комбинация методов 2 + 3 дает достаточную скорость для построения многоугольников. Использование методов 2 + 3 + 1 добавляет проблему потери хороших свойств spобъектов, и моя основная трудность заключается в применении проекций. Я взломал что-то вместе, чтобы спроектировать объект фрейма данных, но он работает довольно медленно. Я думаю, что использование метода 2 + 3 обеспечивает мне достаточное ускорение, пока я не смогу получить изломы от использования метода 2 + 3 + 1.

Каждый должен изучить переход на пакет sf (пространственные объекты) вместо sp. Это значительно быстрее (в данном случае 1/60) и проще в использовании. Вот пример чтения в shp и печати через ggplot2.

Примечание: вам нужно переустановить ggplot2 из самой последней сборки на github (см. Ниже)

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 имеют очень высокое пространственное разрешение береговой линии. Если вам это не нужно, вы можете использовать более обобщенный набор данных. Подходы ialm очень интересны, но простой альтернативой является использование данных 'wrld_simpl', поставляемых с 'maptools'

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