0. 数据及软件包准备
- 随机生成的兴趣点
- 全球大洲尺度shapefile
- 非洲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函数裁剪后的结果图
