在R里面根据shp文件进行点的抠图

[toc]

在R中根据shapefile进行点的抠取

当我们有一个图层文件时候，然后再放些采样点在该图层上，发现，某些采样点落在地图的外面，如下图所示。如果用Arcgis操作，很容易把外围的点抹去，保留图层内的点，那么如果在R里面，实现该操作呢。

image.png

本篇文章，主要介绍在R中实现根据shp文件进行地图点的抠取。

数据准备

首先我们利用广西，云南为案例，随机生成一些散在的点。然后将点的df_point转成SpatialPointsDataFrame格式；同样云南地图也转成SpatialPointsDataFrame格式（我这里是用sf读取，如果用SpatialPointsDataFrame读取，则不用转换），记住这里需要设置地图的投影格式为" proj=longlat ellps=WGS84"

代码语言：javascript复制

# 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对象进行操作。这样一来就可以了。

代码语言：javascript复制

#### 根据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

代码语言：javascript复制

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)

参考

point.in.poly (spatialEco) returns non-intersecting points

arcgis image png point shapefile

0 人点赞