R语言地理空间分析（三）矢量栅格数据操作

sp包矢量数据操作

s1 <- readRDS("Income_schooling.rds")

s2 <- readRDS("dist.rds")

l1 <- readRDS("highway.rds")
#函数bbox获取空间数据边界，sp包
b1<-bbox(s1)
b2<-bbox(s2)
b3<-pmax(b1,b2)
b3[,1]<-pmin(b1[,1],b2[,1])

#aggregate()函数，融合边界,合为一体，raster包
ME<-aggregate(s1)

s1$med<-s1$Income>median(s1$Income)
#以某种属性进行融合
ME.inc<-aggregate(s1,by="med")
#group_by()对s1@data进行以med属性分组，然后summarize()统计各组的中位数
#.$medinc将统计结果赋值给
ME.inc$Income <- s1@data %>% group_by(med) %>%
  summarize(medinc = median(Income)) %>% .$medinc
#计算面积
ME.inc$Area<-gArea(ME.inc,byid=TRUE)/1000000

#intersect()截取相交部分,raster
clp1<-intersect(s1,s2)
#按照输入数据顺序，来决定输出数据类型，clp2为线空间数据
clp2<-intersect(l1,s1)
#clp3为面空间数据
clp3<-intersect(s1,l1)


#union()相交融合，raster包
un1<-union(s1,s2)

sf包矢量数据操作

提取对象几何信息：

st_bbox(): 提取对象的坐标范围
st_area(): 求对象的各要素面积
st_length(): 提取对象的各要素周长
st_geometry(): 提取对象的几何信息
st_drop_geometry(): 删除对象的几何要素，成为数据框
st_coordinates(): 提取各要素几何信息的地理或投影坐标
st_boundary(): 提取各要素的边界，返回线对象
st_centriod(): 提取各要素的质心，返回点对象
st_voronoi(): 生成泰森多边形。
叠置分析：
st_intersection：生成对象内要素两两之间的交集；
st_difference：生成对象内要素两两之间的差集；
st_sym_difference：交集取反操作，即去除对象内两两要素的相交部分，保留不相交部分;
st_crop：用指定地理或投影坐标确定的矩形范围截取矢量对象；
其他函数：
st_simplify：简化对象的拓扑关系；
st_make_grid：在指定地理或投影坐标范围内创建网格；
st_jitter：随机偏移矢量对象各要素的位置；
st_sample：在面数据内部按某种规则选取若干个点（默认随机选取）；
st_line_sample：在线数据上按某种规则选取若干点（默认均匀选取）。

栅格操作

bath2 <- bath * (-1)
bath3 <- bath2 < 0
tm_shape(bath3) + tm_raster(palette = "Greys") + 
  tm_legend(outside = TRUE, text.size = .8)
#重分类reclassify()
#创建一个矩阵，第一列和第二列为重分类的起始值和结束值，第三列为显示值
m <- c(0, 100, 100,  100, 500, 500,  500, 
       1000,  1000, 1000, 11000, 11000)
m <- matrix(m, ncol=3, byrow = T)
#right表示右闭合
bath3 <- reclassify(bath, m, right = T)

#两个栅格操作
elevation<-elev-bath
#平滑处理focal()
f1<-focal(elevation,w=matrix(1,nrow = 11,ncol = 11),fun=mean)

m  <- matrix(c(1,1,1,1,0,1,1,1,1)/8,nrow = 3) 
f2 <- focal(elevation, w=m, fun=sum)
#Sobel过滤器
Sobel <- matrix(c(-1,0,1,-2,0,2,-1,0,1) / 4, nrow=3) 
f3    <- focal(elevation, w=Sobel, fun=sum)

#单元格合并aggregate(),raster
z1<-aggregate(elevation,fact=2,fun=mean,expand=TRUE)
#res()用于查看像素大小
res(z1)