空间多边形数据帧 R 中每个多边形的边界框

Question

如何获取polys中每个多边形的bbox？

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

只给我lonlat，但我想为每个多边形获取lat1lat2和lon1lon2。

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

空间面数据框有几个槽。 @data是数据框，@polygons是面。您可以先尝试str(polys@polygons)看看您得到了什么。如果是多边形列表，则使用sp::bbox函数lapply它

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

Answer 2

之前的答案很好，不幸的是，在大型数据集(数百万个多边形(上对我来说太慢了。

使用 Rcpp 可以大大加快速度(>20 倍(，下面的函数是从 spatDataManagement 包中提取的，您可以将其放入 cpp 文件中并在其上运行sourceCpp以获取该功能，或者您可以安装spatDataManagement并直接使用该函数。此外，它很好地将所有内容放在具有足够列名的 data.frame 中：

> 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

法典：

#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

返回可绘制边界框(以空间多边形数据帧的形式(而不是lapply(polys@polygons, bbox)返回的矩阵列表：

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