[toc]
在R中根据shapefile进行点的抠取
当我们有一个图层文件时候，然后再放些采样点在该图层上，发现，某些采样点落在地图的外面，如下图所示。如果用Arcgis操作，很容易把外围的点抹去，保留图层内的点，那么如果在R里面，实现该操作呢。
image.png
本篇文章，主要介绍在R中实现根据shp文件进行地图点的抠取。
数据准备
首先我们利用广西，云南为案例，随机生成一些散在的点。然后将点的df_point转成SpatialPointsDataFrame格式；同样云南地图也转成SpatialPointsDataFrame格式（我这里是用sf读取，如果用SpatialPointsDataFrame读取，则不用转换），记住这里需要设置地图的投影格式为"+proj=longlat +ellps=WGS84"
# point data
set.seed(124)
df_point=tibble(x=rnorm(100,101,2.3),
y=rnorm(100,24,2.3))
# shapefile
Yun= read_sf("https://geo.datav.aliyun.com/areas_v2/bound/530000.json")%>%
st_transform(., 4326)
# plot
ggplot()+
geom_sf(data=Yun,fill=NA,size=0.2)+
geom_point(data = df_point,aes(x,y))
抠取操作
主要借助于point.in.poly函数对两个SpatialPointsDataFrame对象进行操作。这样一来就可以了。
#### 根据shp文件进行抠图
library(rgdal)
spg = df_point
# 1)point change to SpatialPixelsDataFrame
coordinates(spg) = ~ x + y
proj4string(spg) = CRS("+proj=longlat +ellps=WGS84")
# 2)SHP change to SpatialPixelsDataFrame
Yun_shp = as(Yun, 'Spatial')
proj4string(Yun_shp) = CRS("+proj=longlat +ellps=WGS84")
library(spatialEco)
library(ggspatial)
# 3) intersect points in polygon
df_overlap_sp = point.in.poly(spg, Yun_shp)
# convert to data frame, keeping your data
df_overlap = as.data.frame(df_overlap_sp) %>% na.omit()
#plot
ggplot()+
geom_sf(data=Yun,fill=NA,size=0.2)+
geom_point(data=df_overlap,aes(coords.x1,coords.x2))
# change to sf
df_sf = st_as_sf(df_overlap,coords = c("coords.x1","coords.x2")) %>%
st_set_crs(4326)
ggplot()+
geom_sf(data=Yun,fill=NA,size=0.2)+
geom_sf(data=df_sf,fill=NA,size=0.2)
image.png
image.png
library(ncdf4)
library(rgdal)
library(gdalUtils)
library(raster)
library(rasterVis)
library(sf)
library(exactextractr)
library(tidyverse)
rm(list = ls())
filter=dplyr::filter
select=dplyr::select
set.seed(124)
df_point=tibble(x=rnorm(10,101,2.3),
y=rnorm(10,24,2.3))
# shapefile
Yun= read_sf("https://geo.datav.aliyun.com/areas_v2/bound/530000_full.json")%>%
st_transform(., 4326)
# plot
ggplot()+
geom_sf(data=Yun,fill=NA,size=0.2)+
geom_point(data = df_point,aes(x,y))
#### 根据shp文件进行抠图
library(rgdal)
spg = df_point
# 1)point change to SpatialPixelsDataFrame
coordinates(spg) = ~ x + y
proj4string(spg) = CRS("+proj=longlat +ellps=WGS84")
# 2)SHP change to SpatialPixelsDataFrame
Yun_shp = as(Yun, 'Spatial')
proj4string(Yun_shp) = CRS("+proj=longlat +ellps=WGS84")
library(spatialEco)
library(ggspatial)
# 3) intersect points in polygon
df_overlap_sp = point.in.poly(spg, Yun_shp)
# convert to data frame, keeping your data
df_overlap = as.data.frame(df_overlap_sp) %>% na.omit()
#plot
ggplot()+
geom_sf(data=Yun,fill=NA,size=0.2)+
geom_point(data=df_overlap,aes(coords.x1,coords.x2))
# change to sf
df_sf = st_as_sf(df_overlap,coords = c("coords.x1","coords.x2")) %>%
st_set_crs(4326)
ggplot()+
geom_sf(data=Yun,fill=NA,size=0.2)+
geom_sf(data=df_sf,fill=NA,size=0.2)
df=Yun %>% select(adcode,name,subFeatureIndex,geometry) %>%
mutate(val=ifelse(subFeatureIndex %in% c(df_overlap$subFeatureIndex),1,0))
ggplot()+
geom_sf(data=df,fill=NA,size=0.2)
library(raster)
shape = as(df, 'Spatial')
proj4string(Yun_shp) = CRS("+proj=longlat +ellps=WGS84")
r = raster(shape, res=0.05)
shape_r = rasterize(shape, r, "val")
plot(shape_r)
plot(shape,add=T)
