Boîte englobante pour chaque polygone dans SpatialPolygonsDataFrame R

Question

Comment puis-je obtenir le bbox pour chaque polygone dans polys?

pp <- cbind(coordinates(polys),as.data.frame(polys)) 

me donne lonlat seulement, mais je voudrais obtenir lat1lat2 et lon1lon2 pour chaque polygone.

polys=new("SpatialPolygonsDataFrame" 
     , data = structure(list(NAMRB_EN = structure(c(6L, 45L, 2L, 41L, 31L, 
    3L, 40L, 14L, 42L, 7L, 26L, 12L, 38L, 25L, 36L, 9L, 39L, 27L, 
    32L, 19L, 43L, 21L, 15L, 22L, 20L, 9L, 17L, 11L, 33L, 44L, 37L, 
    13L, 8L, 5L, 18L, 30L, 16L, 10L, 1L, 29L, 34L, 23L, 24L, 28L, 
    4L, 35L), .Label = c("Albany", "Arctic Ocean Seaboard", "Arnaud", 
    "Atlantic Ocean Seaboard", "Attawapiskat", "Back", "Baleine", 
    "Broadback", "Churchill", "Columbia", "Eastmain", "Feuilles", 
    "Fraser", "George", "Grande Baleine", "Harricanaw", "Hayes", 
    "Hudson Bay Seaboard", "Koksoak", "La Grande", "Little Mecatina", 
    "Mackenzie", "Mississippi", "Moose", "Naskaupi", "Nass", "Natashquan", 
    "Nelson", "Nottaway", "Pacific Ocean Seaboard", "Povungnituk", 
    "Romaine", "Rupert", "Saint John", "Saint Lawrence", "Seal", 
    "Severn", "Skeena", "St.-Augustin", "Stikine", "Taku", "Thelon", 
    "Wannock", "Winisk", "Yukon"), class = "factor"), NAODA_EN = structure(c(1L, 
    5L, 1L, 5L, 4L, 4L, 5L, 4L, 4L, 4L, 5L, 4L, 5L, 2L, 4L, 2L, 2L, 
    2L, 2L, 4L, 5L, 2L, 4L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 4L, 
    4L, 4L, 5L, 4L, 5L, 4L, 4L, 2L, 3L, 4L, 4L, 2L, 2L), .Label = c("Arctic Ocean", 
    "Atlantic Ocean", "Gulf of Mexico", "Hudson Bay", "Pacific Ocean" 
    ), class = "factor"), NAMRB_FR = structure(c(4L, 45L, 19L, 41L, 
    31L, 2L, 40L, 12L, 42L, 5L, 27L, 10L, 38L, 26L, 36L, 7L, 39L, 
    28L, 32L, 16L, 43L, 18L, 13L, 23L, 17L, 7L, 15L, 9L, 33L, 44L, 
    37L, 11L, 6L, 3L, 22L, 21L, 14L, 8L, 1L, 30L, 34L, 24L, 25L, 
    29L, 20L, 35L), .Label = c("Albany", "Arnaud", "Attawapiskat", 
    "Back", "Baleine", "Broadback", "Churchill", "Columbia", "Eastmain", 
    "Feuilles", "Fraser", "George", "Grande Baleine", "Harricanaw", 
    "Hayes", "Koksoak", "La Grande", "Little Mecatina", "Littoral de l'océan Arctique", 
    "Littoral de l'océan Atlantique", "Littoral de l'océan Pacifique", 
    "Littoral de la Baie d'Hudson", "Mackenzie", "Mississippi", "Moose", 
    "Naskaupi", "Nass", "Natashquan", "Nelson", "Nottaway", "Povungnituk", 
    "Romaine", "Rupert", "Saint-Jean", "Saint-Laurent", "Seal", "Severn", 
    "Skeena", "St.-Augustin", "Stikine", "Taku", "Thelon", "Wannock", 
    "Winisk", "Yukon"), class = "factor"), NAODA_FR = structure(c(3L, 
    5L, 3L, 5L, 1L, 1L, 5L, 1L, 1L, 1L, 5L, 1L, 5L, 4L, 1L, 4L, 4L, 
    4L, 4L, 1L, 5L, 4L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 5L, 1L, 
    1L, 1L, 5L, 1L, 5L, 1L, 1L, 4L, 2L, 1L, 1L, 4L, 4L), .Label = c("Baie d'Hudson", 
    "Golfe de Mexique", "Océan Arctique", "Océan Atlantique", "Océan Pacifique" 
    ), class = "factor")), .Names = c("NAMRB_EN", "NAODA_EN", "NAMRB_FR", 
    "NAODA_FR"), row.names = 0:45, class = "data.frame") 
     , polygons = list(<S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>, 
     <S4 object of class structure("Polygons", package = "sp")>) 
     , plotOrder = c(3L, 24L, 35L, 46L, 44L, 2L, 42L, 38L, 45L, 36L, 26L, 9L, 32L, 
    1L, 20L, 39L, 27L, 31L, 43L, 25L, 16L, 7L, 30L, 40L, 6L, 15L, 
    34L, 13L, 12L, 41L, 28L, 8L, 23L, 29L, 5L, 10L, 37L, 11L, 14L, 
    33L, 4L, 22L, 18L, 19L, 17L, 21L) 
     , bbox = structure(c(-152.812332679775, 40.3769750107632, -52.6362915039062, 
    83.1106262207029), .Dim = c(2L, 2L), .Dimnames = list(c("x", 
    "y"), c("min", "max"))) 
     , proj4string = new("CRS" 
     , projargs = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0" 
    ) 
    ) 

Answer 1

La zone de données de polygone spatial a quelques emplacements. @data est la trame de données, @polygons est les polygones. Vous pouvez d'abord essayer str(polys@polygons) pour voir ce que vous obtenez. Si elle est une liste de polygones, puis lapply avec sp::bbox fonction

require(sp) 
lapply(polys@polygons, bbox) 

Answer 2

Réponse précédente est bien, malheureusement, il était beaucoup trop lent pour moi sur de grands ensembles de données (en millions de polygones).

L'utilisation CRPP permet d'accélérer beaucoup (> 20 fois), l'au-delà de la fonction est extrait du paquet spatDataManagement, vous pouvez le mettre dans un fichier cpp et exécutez sourceCpp sur elle pour obtenir la fonction ou vous pouvez installer spatDataManagement et utilisez la fonction directement. En outre, il met bien tout en data.frame avec les noms de colonnes adéquates:

> head(GetBBoxes(phillyCensusTracts)) 
     minX  minY  maxX  maxY 
1 -74.97678 40.07721 -74.95576 40.09781 
2 -75.00539 40.09937 -74.96408 40.12390 
3 -75.14641 39.92889 -75.13040 39.96294 
4 -75.00942 40.05475 -74.99043 40.06916 
5 -75.03647 40.05078 -75.02171 40.06731 
6 -74.98463 40.07198 -74.97151 40.09381 

code:

#include <Rcpp.h> 
#include <string> 

using namespace Rcpp; 

//' Get bounding box for each polygon/polyline 
//' @description Gets the bounding box (range of x and y) for each item of a SpatialLines/Polygons (DataFrame or not) object. 
//' It is equivalent to applying the \code{sp::bbox} function to all polygons with lapply but it is 
//' simpler to use and much faster (even on a toy example of a few hundred polygons but designed for datasets with millions, then easily gets > 20x faster). 
//' It is an important function to speed up other more complex comparisons between sp objects such as over(). 
//' @param x A SpatialPolygons[DataFrame] or SpatialLines[DataFrame] 
//' @return A matrix of same number of columns than x and 4 columns : xmin, ymin, xmax, ymax 
//' @export 
//' @examples 
//' data(phillyCensusTract) 
//' ## simple use 
//' system.time(bboxes <- GetBBoxes(phillyCensusTracts)) 
//' #=> 0.002 seconds 
//' ## same using bbox 
//' system.time(bboxesRef <- matrix(unlist(lapply(phillyCensusTracts@polygons,bbox)),ncol=4,byrow=TRUE)) 
//' #=> 0.021 seconds 
// [[Rcpp::export]] 
SEXP GetBBoxes(SEXP x){ 
    // determines object type and adapts the search of coordinates 
    S4 obj(x) ; 
    std::string nameList; 
    std::string nameSubList; 
    if(Rf_inherits(x, "SpatialLines") || Rf_inherits(x, "SpatialLinesDataFrame")){ 
     nameList = "lines"; 
     nameSubList = "Lines"; 
    }else if(Rf_inherits(x, "SpatialPolygons") || Rf_inherits(x, "SpatialPolygonsDataFrame")){ 
     nameList = "polygons"; 
     nameSubList = "Polygons"; 
    }else{ 
     ::Rf_error("In GetBBoxes, class must be Spatial[Polygons|Lines][DataFrame]"); 
    } 
    List a = obj.slot(nameList); 

    // count items 
    int nPol = a.length(); 
    NumericMatrix bboxes(nPol,4); 

    // get the range 
    for(int iPol = 0;iPol < nPol;iPol++){ 
     S4 pol = a(iPol); 
     List b = pol.slot(nameSubList); 

     double minX = std::numeric_limits<double>::infinity(); 
     double maxX = -std::numeric_limits<double>::infinity(); 
     double minY = std::numeric_limits<double>::infinity(); 
     double maxY = -std::numeric_limits<double>::infinity(); 

     for(int iSP = 0; iSP < b.length(); iSP++){ 
      S4 subPol = b(iSP); 
      NumericMatrix coords = subPol.slot("coords"); 
      // X 
      NumericVector rangeX = range(coords(_,0)); 
      if(rangeX(0)<minX) minX = rangeX(0); 
      if(rangeX(1)>maxX) maxX = rangeX(1); 
      // Y 
      NumericVector rangeY = range(coords(_,1)); 
      if(rangeY(0)<minY) minY = rangeY(0); 
      if(rangeY(1)>maxY) maxY = rangeY(1); 
     } 

     bboxes(iPol,0) = minX; 
     bboxes(iPol,1) = minY; 
     bboxes(iPol,2) = maxX; 
     bboxes(iPol,3) = maxY; 
    } 
    Rcpp::DataFrame BBoxes = Rcpp::DataFrame::create(Rcpp::Named("minX")=bboxes(_,0), 
      Rcpp::Named("minY")=bboxes(_,1), 
      Rcpp::Named("maxX")=bboxes(_,2), 
      Rcpp::Named("maxY")=bboxes(_,3)); 

    return BBoxes; 
} 

Answer 3

Pour revenir boîtes englobantes traçable (sous la forme d'un SpatialPolygonsDataFrame) au lieu d'une liste de matrices lapply(polys@polygons, bbox) retours:

spatial_bboxes <- function(polygons) { 
    individual_bb <- function(polygon, projection) { 
    polygon <- sp::SpatialPolygons(list(polygon), proj4string = projection) 
    spatial_bbox <- as(raster::extent(polygon), "SpatialPolygons") 
    spatial_bbox <- spatial_bbox@polygons[[1]] 
    spatial_bbox@ID <- polygon@polygons[[1]]@ID 
    return(spatial_bbox) 
    } 
    polys <- lapply(polygons@polygons, individual_bb, polygons@proj4string) 
    spatial_polys <- sp::SpatialPolygons(polys, proj4string = polygons@proj4string) 
    spatial_polys_df <- sp::SpatialPolygonsDataFrame(spatial_polys, polygons@data) 
    return(spatial_polys_df) 
}