2017-04-06 191 views
1

如何獲得polys中每個多邊形的bboxSpatialPolygonsDataFrame中每個多邊形的邊界框R

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

給我只能lonlat,但我想獲得lat1lat2lon1lon2爲每個多邊形。

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" 
    ) 
    ) 
+0

我不能讓code123的數據加載,但我的一個例子,提出了@代碼李哲源ZheyuanLi工作的罰款。 – G5W

+0

@李哲源ZheyuanLi真棒。請提供以供我接受。 – code123

回答

2

空間多邊形數據幀有幾個插槽。 @data是數據幀,@polygons是多邊形。你可以先嚐試str([email protected])看看你得到了什麼。如果是多邊形的列表,然後lapplysp::bbox功能

require(sp) 
lapply([email protected], bbox) 
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([email protected],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; 
} 
0

要返回plottable邊框(在SpatialPolygonsDataFrame的形式),而不是名單這lapply([email protected], 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 <- [email protected][[1]] 
    [email protected] <- [email protected][[1]]@ID 
    return(spatial_bbox) 
    } 
    polys <- lapply([email protected], individual_bb, [email protected]) 
    spatial_polys <- sp::SpatialPolygons(polys, proj4string = [email protected]) 
    spatial_polys_df <- sp::SpatialPolygonsDataFrame(spatial_polys, [email protected]) 
    return(spatial_polys_df) 
}