如何在 SPDF 中使用距离对多边形进行子集化
How to subset polygons using distance in a SPDF
好的。所以，我正在学习并提出了一个非常相似的问题，当时我还没有学到什么 Spatial Polygons Data Frame was about 10 days ago: .
现在，我发现了 SPDF 和等值线图的魔力，并且在本质上遇到了相同的问题，但文件类型不同。我仍在思考 S4 对象，我无法弄清楚如何从我的数据集中对某些迷你多边形 MUNICIPI 进行子集化。
进入正题！
上下文：
我有一个包含这些数据的 SPDF：
目标：
我想对距离海岸一定距离内的所有 MUNICIPI 进行子集化。假设@urbandatascientist 在他对我的第一个问题的回答中设置了 20 公里，并创建了一个等值线图，例如 upper_trees 的子集 MUNICIPI.
从 post 我还有 list.of.polygon.boundaries 我们将从中减去 MUNICIPI 坐标。
一旦我对沿海地区 MUNICIPI 进行了子集化，我希望地图看起来像 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(legend.outside =TRUE)
概述
类似于Select raster in ggplot near coastline的答案，解决方案涉及以下步骤：
正在计算地中海西部 shape file 的巴利阿里（伊比利亚海）和西部盆地部分的边界坐标。
计算 MUNICIPI 中每个多边形的质心，从 OP 的 link 到她的 Google Drive 文件夹，其中包含形状文件的 .zip 文件.
计算前两点与子集MUNICIPI之间的距离，显示其质心到西部盆地的距离小于或等于20公里的多边形。
筛选西部盆地坐标
我没有计算每个质心到 western.basin.polygon.coordinates 内每个坐标对的距离，而是只包括 western.basin.polygon.coordinates 内的坐标，其纬度点在加泰罗尼亚东海岸之间（含）。
作为参考，我使用Peniscola, Catalonia and Cerbere, France的纬度点。通过仅保留位于加泰罗尼亚东海岸的坐标对，western.basin.polygon.coordinates 和 MUNICIPI 中每个质心之间的距离计算在大约 4 分钟内完成。
然后我将 MUNICIPI 内的多边形索引存储在 less.than.or.equal.to.max.km 中，其质心距离小于或等于 20 公里 - TRUE/FALSE 值的逻辑向量。使用 leaflet 地图，我展示了如何对 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
好的。所以，我正在学习并提出了一个非常相似的问题，当时我还没有学到什么 Spatial Polygons Data Frame was about 10 days ago:
现在，我发现了 SPDF 和等值线图的魔力，并且在本质上遇到了相同的问题，但文件类型不同。我仍在思考 S4 对象，我无法弄清楚如何从我的数据集中对某些迷你多边形 MUNICIPI 进行子集化。
进入正题！
上下文：
我有一个包含这些数据的 SPDF：
目标：
我想对距离海岸一定距离内的所有 MUNICIPI 进行子集化。假设@urbandatascientist 在他对我的第一个问题的回答中设置了 20 公里，并创建了一个等值线图，例如 upper_trees 的子集 MUNICIPI.
从 list.of.polygon.boundaries 我们将从中减去 MUNICIPI 坐标。
一旦我对沿海地区 MUNICIPI 进行了子集化，我希望地图看起来像 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(legend.outside =TRUE)
概述
类似于Select raster in ggplot near coastline的答案，解决方案涉及以下步骤：
正在计算地中海西部 shape file 的巴利阿里（伊比利亚海）和西部盆地部分的边界坐标。
计算
MUNICIPI中每个多边形的质心，从 OP 的 link 到她的 Google Drive 文件夹，其中包含形状文件的 .zip 文件.计算前两点与子集
MUNICIPI之间的距离，显示其质心到西部盆地的距离小于或等于20公里的多边形。
筛选西部盆地坐标
我没有计算每个质心到 western.basin.polygon.coordinates 内每个坐标对的距离，而是只包括 western.basin.polygon.coordinates 内的坐标，其纬度点在加泰罗尼亚东海岸之间（含）。
作为参考，我使用Peniscola, Catalonia and Cerbere, France的纬度点。通过仅保留位于加泰罗尼亚东海岸的坐标对，western.basin.polygon.coordinates 和 MUNICIPI 中每个质心之间的距离计算在大约 4 分钟内完成。
然后我将 MUNICIPI 内的多边形索引存储在 less.than.or.equal.to.max.km 中，其质心距离小于或等于 20 公里 - TRUE/FALSE 值的逻辑向量。使用 leaflet 地图，我展示了如何对 MUNICIPI 进行子集化以仅可视化那些在 less.than.or.equal.to.max.km.
# 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