[toc]

在R中根據(jù)shapefile進(jìn)行點(diǎn)的摳取

當(dāng)我們有一個(gè)圖層文件時(shí)候，然后再放些采樣點(diǎn)在該圖層上，發(fā)現(xiàn)，某些采樣點(diǎn)落在地圖的外面，如下圖所示。如果用Arcgis操作，很容易把外圍的點(diǎn)抹去，保留圖層內(nèi)的點(diǎn)，那么如果在R里面，實(shí)現(xiàn)該操作呢。

image.png

本篇文章，主要介紹在R中實(shí)現(xiàn)根據(jù)shp文件進(jìn)行地圖點(diǎn)的摳取。

數(shù)據(jù)準(zhǔn)備

首先我們利用廣西，云南為案例，隨機(jī)生成一些散在的點(diǎn)。然后將點(diǎn)的df_point轉(zhuǎn)成SpatialPointsDataFrame格式；同樣云南地圖也轉(zhuǎn)成SpatialPointsDataFrame格式（我這里是用sf讀取，如果用SpatialPointsDataFrame讀取，則不用轉(zhuǎn)換），記住這里需要設(shè)置地圖的投影格式為"+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函數(shù)對(duì)兩個(gè)SpatialPointsDataFrame對(duì)象進(jìn)行操作。這樣一來(lái)就可以了。

#### 根據(jù)shp文件進(jìn)行摳圖
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))


#### 根據(jù)shp文件進(jìn)行摳圖
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