r - 将地理引用数据与 R 中的形状文件匹配

Question

我有一个地理引用事件数据集，格式为:

LONGITUDE LATITUDE VAR1
33.4      4.4      5
33.4      4.4      3
33.4      4.4      1
30.4      4.2      2
28.4      5.1      2

它计算地理引用事件中的死亡人数。除此之外，我有一个国家省份的形状文件，如下所示:

> str(shapefile)
'data.frame':   216 obs. of  5 variables:
 $ CONSTI_COD: num  1 2 3 4 5 6 7 8 9 10 ...
 $ Area      : num  20 11.7 10.7 223.3 38.7 ...
 $ PROVINCE_NAME  : Factor w/ 216 levels "CENTRAL","COAST",..: 4 4 4 4 4 4 4 4 2 2 ...
 $ Shape_Leng: num  0.193 0.152 0.201 0.872 0.441 ...
 $ Shape_Area: num  0.001628 0.000947 0.000867 0.018135 0.003145 ...

 ..@ polygons   :List of 216
   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
   .. .. .. ..@ Polygons :List of 1
   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
   .. .. .. .. .. .. ..@ labpt  : num [1:2] 36.9 -1.3
   .. .. .. .. .. .. ..@ area   : num 0.00163
   .. .. .. .. .. .. ..@ hole   : logi FALSE
   .. .. .. .. .. .. ..@ ringDir: int 1
   .. .. .. .. .. .. ..@ coords : num [1:151, 1:2] 36.8 36.8 36.8 36.9 36.9 ...
   .. .. .. ..@ plotOrder: int 1
   .. .. .. ..@ labpt    : num [1:2] 36.9 -1.3
   .. .. .. ..@ ID       : chr "0"
   .. .. .. ..@ area     : num 0.00163
 [...etc]

我需要做的是将事件数据放在省份内，即根据坐标将第四列添加到第一个数据框中，该列说明每个事件发生在哪个省份。所以我会有这样的事情:

LONGITUDE LATITUDE VAR1 PROVINCE
33.4      4.4      5    CENTRAL
33.4      4.4      3    CENTRAL
33.4      4.4      1    CENTRAL
30.4      4.2      2    COAST
28.4      5.1      2    COAST

这可能吗？我想我前段时间找到了一篇解释如何做到这一点的帖子(在 Stack Overflow 之外)，但我现在找不到了。

谢谢!

(对不起，如果这里有类似的问题。我进行了搜索，但没有找到答案，也许是因为我真的不知道我在找什么。我真的很感激类似帖子的链接.)

Answer 1

您所谈论的是“空间连接”(或“空间交集”或“叠加”)。在 over 的帮助下，这非常简单。来自 sp 的函数包裹。

这是一个例子。

首先，让我们下载并导入世界各国的多边形 shapefile。

download.file(paste0('http://www.naturalearthdata.com/http//',
                     'www.naturalearthdata.com/download/110m/cultural/',
                     'ne_110m_admin_0_countries.zip'), 
              f <- tempfile())
unzip(f, exdir=tempdir())
library(rgdal)
countries <- readOGR(tempdir(), 'ne_110m_admin_0_countries')

现在我们将创建一些落在多边形 shapefile 范围内的随机坐标数据。然后我们定义列 x和 y如 coordinates ，并分配与多边形相同的 CRS(尽管您的数据可能不是这种情况；请确保分配正确的坐标系)。

pts <- data.frame(x=runif(10, -180, 180), y=runif(10, -90, 90),
                  VAR1=LETTERS[1:10])
coordinates(pts) <- ~x+y  # pts needs to be a data.frame for this to work
proj4string(pts) <- proj4string(countries)

plot(countries)
points(pts, pch=20, col='red')

现在我们可以执行空间叠加:

over(pts, countries)$admin

#  [1] <NA>      <NA>      Turkey    <NA>     
#  [5] Macedonia <NA>      China     Argentina
#  [9] <NA>      Canada   
# 177 Levels: Afghanistan Albania ... Zimbabwe

请注意，在这种情况下，一些随机点落在海洋中(即多边形外)。当与多边形对象相交时，这些点返回 NA。

现在我们 cbind pts 所需的属性:

cbind.data.frame(pts, country=over(pts, countries)$admin)

#             x          y VAR1   country
# 1   -52.59404 -37.422879    A      <NA>
# 2   -33.88867 -40.194482    B      <NA>
# 3    38.84383  37.272460    C    Turkey
# 4   -84.04949   7.118878    D      <NA>
# 5    20.98272  40.920470    E Macedonia
# 6  -155.32951 -37.612497    F      <NA>
# 7    99.40166  38.630049    G     China
# 8   -61.84025 -27.412885    H Argentina
# 9   -37.65287  -3.666080    I      <NA>
# 10 -112.81197  59.959475    J    Canada