r - 如何在 SPDF 中使用距离对多边形进行子集化

Question

好的。所以，我正在学习并提出了一个非常相似的问题，当时我还没有学到什么是 Spatial Polygons Data Frame。大约 10 天前:Select raster in ggplot near coastline .

现在，我发现了 SPDF 和 choropleth map 的魔力，本质上是相同的问题，但文件类型不同。我仍然在思考 S4 对象，我无法弄清楚如何从我的数据集中对某些迷你多边形 MUNICIPI 进行子集化。

进入正题!

上下文: 我有一个包含这些数据的 SPDF:

目标:

我想对距离海岸一定距离内的所有 MUNICIPI 进行子集化。假设@urbandatascientist 在他对我的第一个问题的回答中设置了 20 公里，并创建了一个等值线图，例如，upper_trees 和子集 MUNICIPI。

来自Select raster in ggplot near coastline我还有 list.of.polygon.boundaries，我们将从中减去 MUNICIPI 坐标。

Data is here .

一旦我对沿海 MUNICIPI 进行了子集化，我希望 map 看起来会是 something like the green shaded region here.我还尝试确保 list.of.polygon.boundaries 之间的坐标相同。

任何线索或想法将不胜感激!

到目前为止，这是我使用 upper_trees 作为示例的整个区域的叶绿体图:

tm_shape(catasense2)+ tm_fill(col="upper_trees",n=8,style="quantile")+ tm_layout(图例.outside =TRUE)

Answer 1

概览

类似于 Select raster in ggplot near coastline 的答案，解决方案包括以下步骤:

• 计算巴利阿里(伊比利亚海)和西部盆地部分边界的坐标 shape file地中海西部。

• 根据 OP 到她的 Google 云端硬盘文件夹的链接计算 MUNICIPI 中每个多边形的质心，该文件夹包含形状文件的 .zip 文件。

• 计算前两点与子集MUNICIPI之间的距离，以显示其质心到西部盆地的距离小于或等于20公里的多边形。

过滤西部盆地的坐标

我没有计算每个质心到 western.basin.polygon.coordinates 中每个坐标对的距离，我只在 western.basin.polygon.coordinates 中包含了坐标其纬度点在加泰罗尼亚东海岸之间(含)。

作为引用，我使用 Peniscola, Catalonia 的纬度点和 Cerbere, France .通过仅保留位于加泰罗尼亚东海岸的坐标对，western.basin.polygon.coordinates 与 MUNICIPI 中每个质心之间的距离计算大约在 ~4 内完成分钟。

然后我将多边形的索引存储在 MUNICIPI 中，其质心距离小于或等于 20 公里 less.than.or.equal.to.max.km - TRUE/FALSE 值的逻辑向量。使用 leaflet map ，我展示了如何对 MUNICIPI 进行子集化以仅可视化那些在 less.than.or.equal.to.max.km 内包含 TRUE 值的多边形。

# load necessary packages
library( geosphere )
library( leaflet )
library( sf )

# load necessary data
unzip( zipfile = "catasense2-20180327T023956Z-001.zip" )

# transfrom zip file contents
# into a sf data frame
MUNICIPI <-
  read_sf(
    dsn = paste0( getwd(), "/catasense2" )
    , layer = "catasense2"
  )


# map data values to colors
my.color.factor <- 
  colorFactor( palette = "viridis", domain = MUNICIPI$uppr_tr )

# Visualize
leaflet( data = MUNICIPI ) %>%
  setView( lng = 1.514619
           , lat = 41.875227
           , zoom = 8 ) %>%
  addTiles() %>%
  addPolygons( color = ~my.color.factor( x = uppr_tr ) )

# unzip the western basin zip file
# unzip( zipfile = "iho.zip" )

# create sf data frame
# of the western basin
western.basin <-
  read_sf( dsn = getwd()
           , layer = "iho"
           , stringsAsFactors = FALSE )

# filter the western.basin
# to only include those bodies of water
# nearest catalonia
water.near.catalonia.condition <-
  which( 
    western.basin$name %in% 
      c( "Balearic (Iberian Sea)"
         , "Mediterranean Sea - Western Basin" )
    )

western.basin <-
  western.basin[ water.near.catalonia.condition, ]

# identify the coordinates for each of the
# remaining polygons in the western.basin
# get the boundaries of each
# polygon within the western basin
# and retain only the lon (X) and lat (Y) values
western.basin.polygon.coordinates <-
  lapply( 
    X = western.basin$geometry
    , FUN = function( i )
      st_coordinates( i )[ , c( "X", "Y" ) ]
  )

# find the centroid
# of each polygon in MUNICIPI
MUNICIPI$centroids <-
  st_centroid( x = MUNICIPI$geometry )

# Warning message:
#   In st_centroid.sfc(x = MUNICIPI$geometry) :
#   st_centroid does not give correct centroids for longitude/latitude data

# store the latitude
# of the western (Peniscola, Catalonia) and eastern (Cerbere, France) 
# most parts of the eastern coast of Catalonia
east.west.lat.range <- 
  setNames( object = c( 40.3772185, 42.4469981 )
            , nm = c( "east", "west" )
  )

# set the maximum distance
# allowed between a point in df
# and the sea to 20 kilometers
max.km <- 20

# store the indices in the
# western basin polygon coordinates whose
# lat points
# fall in between (inclusive) 
# of the east.west.lat.range
satisfy.east.west.lat.max.condition <-
  lapply(
    X = western.basin.polygon.coordinates
    , FUN = function(i)
      which(
        i[, "Y"] >= east.west.lat.range["east"] &
          i[, "Y"] <= east.west.lat.range["west"]
      )
  )

# filter the western basin polygon coordinate
# by the condition
western.basin.polygon.coordinates <-
  mapply(
    FUN = function( i, j )
      i[ j, ]
    , western.basin.polygon.coordinates
    , satisfy.east.west.lat.max.condition
  )

# calculate the distance of each centroid
# in MUNICIPI
# to each point in both western.basin
# polygon boundary coordinates
# Takes ~4 minutes
distance <-
  apply(
    X = do.call( what = "rbind", args = MUNICIPI$centroids )
    , MARGIN = 1
    , FUN = function( i )
      lapply(
        X = western.basin.polygon.coordinates
        , FUN = function( j )
          distGeo(
            p1 = i
            , p2 = j
          ) / 1000 # to transform results into kilometers
      )
  )

# 1.08 GB list is returned
object.size( x = distance)
# 1088805736 bytes

# find the minimum distance value
# for each list in distance
distance.min <-
  lapply(
    X = distance
    , FUN = function( i )
      lapply(
        X = i
        , FUN = function( j )
          min( j )
      )
  )

# identify which points in df
# are less than or equal to max.km
less.than.or.equal.to.max.km <-
  sapply(
    X = distance.min
    , FUN = function( i )
      sapply(
        X = i
        , FUN = function( j )
          j <= max.km
      )
  )

# convert matrix results into
# vector of TRUE/FALSE indices
less.than.or.equal.to.max.km <-
  apply(
    X = less.than.or.equal.to.max.km
    , MARGIN = 2
    , FUN = any
  )

# now only color data that meets
# the less.than.or.equal.to.max.km condition
leaflet( data = MUNICIPI[ less.than.or.equal.to.max.km, ] ) %>%
  setView( lng = 1.514619
           , lat = 41.875227
           , zoom = 8 ) %>%
  addTiles() %>%
  addPolygons( color = ~my.color.factor( x = uppr_tr ) ) 

# end of script #

session 信息

R version 3.4.4 (2018-03-15)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.2

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] sf_0.6-0           leaflet_1.1.0.9000 geosphere_1.5-7   

loaded via a namespace (and not attached):
 [1] mclust_5.4         Rcpp_0.12.16       mvtnorm_1.0-7     
 [4] lattice_0.20-35    class_7.3-14       rprojroot_1.3-2   
 [7] digest_0.6.15      mime_0.5           R6_2.2.2          
[10] plyr_1.8.4         backports_1.1.2    stats4_3.4.4      
[13] evaluate_0.10.1    e1071_1.6-8        ggplot2_2.2.1     
[16] pillar_1.2.1       rlang_0.2.0        lazyeval_0.2.1    
[19] diptest_0.75-7     whisker_0.3-2      kernlab_0.9-25    
[22] R.utils_2.6.0      R.oo_1.21.0        rmarkdown_1.9     
[25] udunits2_0.13      stringr_1.3.0      htmlwidgets_1.0   
[28] munsell_0.4.3      shiny_1.0.5        compiler_3.4.4    
[31] httpuv_1.3.6.2     htmltools_0.3.6    nnet_7.3-12       
[34] tibble_1.4.2       gridExtra_2.3      dendextend_1.7.0  
[37] viridisLite_0.3.0  MASS_7.3-49        R.methodsS3_1.7.1 
[40] grid_3.4.4         jsonlite_1.5       xtable_1.8-2      
[43] gtable_0.2.0       DBI_0.8            magrittr_1.5      
[46] units_0.5-1        scales_0.5.0       stringi_1.1.7     
[49] reshape2_1.4.3     viridis_0.5.0      flexmix_2.3-14    
[52] sp_1.2-7           robustbase_0.92-8  squash_1.0.8      
[55] RColorBrewer_1.1-2 tools_3.4.4        fpc_2.1-11        
[58] trimcluster_0.1-2  DEoptimR_1.0-8     crosstalk_1.0.0   
[61] yaml_2.1.18        colorspace_1.3-2   cluster_2.0.6     
[64] prabclus_2.2-6     classInt_0.1-24    knitr_1.20        
[67] modeltools_0.2-21