R-如何在R中用shapefile掩膜兴趣点

0. 数据及软件包准备

1. 随机生成的兴趣点
2. 全球大洲尺度shapefile
3. 非洲shapefile

# 0. import required packages
require(ggplot2)
require(raster)
# 1. import world and asia shp
setwd('...文件路径...')
world = shapefile('world_continent2.shp')
africa = shapefile('africa.shp')
print(class(world))
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
print(class(africa))
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
#2. Generate Random Points
long = runif(500,min = -180,max = 180)
lat = runif(500,min = -90, max = 83)
df = data.frame(long, lat)

1. 目标

利用Africa shapefile 掩膜裁剪随机生成的世界范围的兴趣点

# 3. 图1 绘图代码
 a = ggplot()+
  geom_polygon(data = world,
               aes(x = long,
                   y = lat,
                   group = group),
               color ='black',
               fill = '#CFCFCF',
               size = 1)+
  geom_polygon(data = africa,
               aes(x = long,
                   y = lat,
                   group = group),
               color = 'black',
               fill = '#CD5C5C',
               size =1)+
  coord_fixed(1.3)+
  geom_point(data = df, aes(x = long,
                            y = lat),
             shape = 21, color = '#1E90FF',
             size = 3)+
  
  geom_point(data = df, aes(x = long,
                            y = lat),
             shape = 21, color = '#1E90FF',
             size = 2.8)+
  
  geom_point(data = df, aes(x = long,
                            y = lat),
             shape = 21, color = '#1E90FF',
             size = 3.2)+
  theme(panel.background = element_rect(fill='white',color = 'black',linetype = 'solid',size=1))+
  xlab('Longitude')+
  ylab('Latitude')
a

图1 随机生成全球范围内的随机点

2. 问题分析

shapefile 格式为SpatialPolygonDataFrame，点坐标为data.frame，传统解决思路为：

方案一

将data.frame 转为 SpatialPoints， 然而利用raster包中的mask 函数进行掩膜

# 4. 将data.frame 转化为 SpatialPoints
coordinates(df) = ~ long+lat
proj4string(df) = crs(world)
df
class       : SpatialPoints 
features    : 500 
extent      : -179.4668, 179.8964, -89.98775, 82.8965  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

# 5. 利用raster 包中的mask 工具裁剪
df_m = mask(df, asia) 
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘mask’ for signature ‘"SpatialPoints", "SpatialPolygonsDataFrame"’

方案一报错

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘mask’ for signature ‘"SpatialPoints", "SpatialPolygonsDataFrame"’

方案二

基于africa.shp 的Extent，采用raster包中的crop函数对兴趣点进行裁剪(图2)。根据图2，我们可以发现基于Extent的裁剪会保留边界外的点。

# 6. 基于Africa Extent 裁剪
df_m = crop(df, extent(africa))
df_m
class       : SpatialPoints 
features    : 216 
extent      : -156.5301, 178.1786, -22.17424, 55.17034  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

df_m.df = as.data.frame(df_m)

#7. Visualize the result
b = ggplot()+
  geom_polygon(data = africa,
               aes(x = long,
                   y = lat,
                   group = group),
               color = 'black',
               fill = '#CD5C5C',
               size = 1)+
  geom_point(data = df_m.df,
             aes(x = long,y = lat),
             shape = 21,
             color = '#1E90FF',
             size = 3)+
  coord_fixed(1.3)+
xlab('Longitude')+
ylab('Latitude')
b

图2 基于Extent 裁剪后的结果图

3. 解决方案-实测有效

基于over函数对兴趣点进行掩膜裁剪。

# 8. 基于over函数对兴趣点进行掩膜裁剪
df_m = over(df, africa)
df_index = which(!is.na(df_m))
df_m = df[df_index,]
df_m.df = as.data.frame(df_m,xy = T)

#9. Visualize the result.
c = ggplot()+
  geom_polygon(data = africa,
               aes(x = long,
                   y = lat,
                   group = group),
               color = 'black',
               fill = '#CD5C5C',
               size = 1)+
  geom_point(data = df_m.df,
             aes(x = long,y = lat),
             shape = 21,
             color = '#1E90FF',
             size = 3)+
  coord_fixed(1.3)+
  xlab('Longitude')+
  ylab('Latitude')

c

图3 基于Over函数裁剪后的结果图