Comment accélérer le tracé des polygones dans R?

Je veux tracer les frontières des pays d'Amérique du Nord sur une image raster représentant une variable, puis superposer les contours au-dessus du tracé à l'aide de R. J'ai réussi à le faire en utilisant des graphiques de base et un treillis, mais il semble que le processus de traçage soit beaucoup trop lent! Je ne l'ai pas encore fait dans ggplot2, mais je doute qu'il s'en tirera mieux en termes de vitesse.

J'ai les données dans un fichier netcdf créé à partir d'un fichier grib. Pour l'instant, j'ai téléchargé les frontières des pays pour le Canada, les États-Unis et le Mexique, qui étaient disponibles dans les fichiers RData de GADM qui lisent dans R en tant qu'objets SpatialPolygonsDataFrame.

Voici du 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

Existe-t-il un moyen d'accélérer le tracé des polygones? Sur le système sur lequel je travaille, le tracé prend plusieurs minutes. Je veux éventuellement créer une fonction qui générera facilement un certain nombre de ces tracés pour inspection, et je pense que je tracerai plusieurs de ces cartes, donc je veux augmenter la vitesse des tracés!

Merci!

J'ai trouvé 3 façons d'augmenter la vitesse de traçage des frontières du pays à partir de fichiers de formes pour R. J'ai trouvé de l'inspiration et du code ici et ici .

(1) Nous pouvons extraire les coordonnées des fichiers de formes pour obtenir la longitude et les latitudes des polygones. Ensuite, nous pouvons les mettre dans un bloc de données avec la première colonne contenant les longitudes et la deuxième colonne contenant les latitudes. Les différentes formes sont séparées par des NA.

(2) Nous pouvons supprimer certains polygones de notre fichier de formes. Le fichier de forme est très, très détaillé, mais certaines formes sont de minuscules îles sans importance (pour mes tracés, en tout cas). Nous pouvons définir un seuil minimal de zone de polygone pour conserver les plus grands polygones.

(3) Nous pouvons simplifier la géométrie de nos formes en utilisant l' algorithme Douglas-Peuker . Les bords de nos formes polygonales peuvent être simplifiés, car ils sont très complexes dans le fichier d'origine. Heureusement, il existe un package,, rgeosqui implémente cela.

Installer:

# 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éthode 1: extraire les coordonnées des fichiers de formes dans un bloc de données et des lignes de tracé

L'inconvénient majeur est que nous perdons ici des informations par rapport à la conservation de l'objet en tant qu'objet SpatialPolygonsDataFrame, comme la projection. Cependant, nous pouvons le transformer en objet sp et ajouter les informations de projection, et c'est toujours plus rapide que de tracer les données d'origine.

Notez que ce code s'exécute très lentement sur le fichier d'origine car il y a beaucoup de formes et le bloc de données résultant fait environ 2 millions de lignes.

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

Méthode 2: supprimer les petits polygones

Il existe de nombreuses petites îles qui ne sont pas très importantes. Si vous vérifiez certains des quantiles des zones pour les polygones, nous voyons que beaucoup d'entre eux sont minuscules. Pour l'intrigue du Canada, je suis passé du tracé de plus de mille polygones à seulement des centaines de polygones.

Quantiles pour la taille des polygones pour le Canada:

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

Méthode 3: simplifiez la géométrie des formes polygonales

Nous pouvons réduire le nombre de sommets dans nos formes polygonales en utilisant la gSimplifyfonction du rgeospackage

Code:

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

Quelques repères:

J'ai utilisé le temps écoulé system.timepour comparer mes temps de traçage. Notez que ce ne sont que les moments pour tracer les pays, sans les courbes de niveau et autres choses supplémentaires. Pour les objets sp, je viens d'utiliser la plotfonction. Pour les objets du bloc de données, j'ai utilisé la plotfonction avec type='l'etlines fonction.

Tracer les polygones originaux du Canada, des États-Unis et du Mexique:

73.009 secondes

En utilisant la méthode 1:

2,449 secondes

En utilisant la méthode 2:

17,660 secondes

En utilisant la méthode 3:

16,695 secondes

En utilisant la méthode 2 + 1:

1,729 secondes

En utilisant la méthode 2 + 3:

0,445 seconde

En utilisant la méthode 2 + 3 + 1:

0.172 secondes

D'autres remarques:

Il semble que la combinaison des méthodes 2 + 3 accélère suffisamment le tracé des polygones. L'utilisation des méthodes 2 + 3 + 1 ajoute le problème de perdre les belles propriétés des spobjets, et ma principale difficulté est d'appliquer des projections. J'ai piraté quelque chose ensemble pour projeter un objet de bloc de données, mais cela fonctionne plutôt lentement. Je pense que l'utilisation de la méthode 2 + 3 me fournit des accélérations suffisantes jusqu'à ce que je puisse me débarrasser de l'utilisation de la méthode 2 + 3 + 1.

Tout le monde devrait envisager de transférer vers le package sf (caractéristiques spatiales) au lieu de sp. Il est nettement plus rapide (1 / 60ème dans ce cas) et plus simple d'utilisation. Voici un exemple de lecture dans un shp et de traçage via ggplot2.

Remarque: vous devez réinstaller ggplot2 à partir de la version la plus récente de github (voir ci-dessous)

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 

Les données GADM ont une résolution spatiale très élevée des côtes. Si vous n'en avez pas besoin, vous pouvez utiliser un ensemble de données plus généralisé. Les approches d'ialm sont très intéressantes, mais une alternative simple consiste à utiliser les données 'wrld_simpl' fournies avec 'maptools'

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