R语言——空间数据分析

大白成长记Jan

已于 2023-10-08 08:38:58 修改

阅读量661

点赞数

分类专栏： R 文章标签：数据分析 r语言数据可视化

于 2023-09-27 22:13:34 首次发布

本文链接：https://blog.csdn.net/weixin_45443016/article/details/133338546

版权

R 专栏收录该内容

17 篇文章 2 订阅

订阅专栏

library(ggplot2)
library(plyr)
library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(sf)
library(ggspatial)
library(tmap)
library(tidyverse)
library(patchwork)
getwd()#确认没有中文字符
setwd(dir="D:\\new D 盘文件\\研究生\\研一下学期课程\\R\\2023.5.18课件\\SpacialAnalysis")

##储存工作目录
wdir <- getwd()
wdir

一、空间数据分析-读取与显示

1.1 矢量文件

1.1.1 shp格式

#.shp格式的一般读取
#install.packages("rgdal")
library(rgdal)
#这段代码是用来查询当前计算机上已安装的支持的矢量数据格式驱动程序列表的。ogrDrivers()是rgdal包中的一个函数，它返回所有可用的OGR驱动程序的元数据，并列出了它们可以打开的文件格式和它们支持的写入格式。
type <- ogrDrivers()
#OGR support is provided by the sf and terra packages among others 
#这段代码将两个字符串wdir和"/data/china"连接，sep参数指之间使用的分隔符,完整路径
dsn <- paste(wdir,"/data/china", sep="")
#ogrListLayers()函数将返回一个包含该数据源中所有图层名称的字符向量。
ogrListLayers(dsn)  
#获取指定数据源中指定图层的元数据信息的。dsn参数是用于指定数据源文件路径的字符串，layer参数是指定data/china中要查询的图层名称,返回元数据
ogrInfo(dsn=dsn, layer="bou2_4p")  
ogrInfo(dsn=dsn, layer="bou2_4l") 
#****************************************#
map_rl <- rgdal::readOGR (dsn=dsn, layer="bou2_4p") #读取.shp文件返回一个
# 写法等同于
#library(rgdal)
#map_rl <- readOGR(dsn = "/path/to/shapefile", layer = "layer_name")
names(map_rl)#map_rl 对象中所有成员的名称
slotNames(map_rl)#map_rl 所有的 S4 成员名称（包括继承的成员）

str(slot(map_rl, "data"))#将成员变量变为str
#S4 是一种面向对象的编程模型，用于创建复杂的数据结构和类。S4 类由多个成员变量（也称为“属性”、“字段”或“特性”）组成，每个成员可以保存不同类型的数据。成员变量通常被称为“Slots”
#在 SpatialPointsDataFrame 类型中，data 就是一个成员变量，用于存储每个点的属性信息。其他常见的 S4 类包括 data.frame、matrix、list 等。
summary(map_rl)
plot(map_rl)
plot(map_rl, axes=T)
#axes是边框
plot(map_rl, axes=T, border="blue", col="gray80")

在这里插入图片描述

1.1.2 json格式

par(mfrow=c(1,1))
library(rgdal)
file_json <- paste(wdir,"/data/", "China.json", sep="")
#file_json是文件路径，encoding="UTF-8"指定文件编码为UTF-8，use_iconv=T表示在读取文件前先转换字符编码。
map_json <- readOGR(file_json,encoding="UTF-8",use_iconv=T)
plot(map_json, axes = TRUE)

在这里插入图片描述

1.2 栅格数据

1.2.1 tiff格式

library(rgdal)
setwd(paste(wdir,"/data/GlobeData", sep=""))

map_tif <- readGDAL("wsiearth.tif")
#summary(tif/shp)返回投影范围与坐标系、data的统计信息等
summary(map_tif)

#地理空间数据文件的元数据包括信息文件格式、投影信息、像素尺寸、波段数等等
GDALinfo("wsiearth.tif")
#****************************************#
ds <- GDAL.open("wsiearth.tif")
#GDAL.open() 返回一个数据集对象，不能直接读，而 readGDAL() 直接返回数据的矩阵或数组形式

#获取一个数据集对象 ds 中的颜色表（color table）
coltable <- getColorTable(ds)

#显示图像，xlim 和 ylim：指定坐标范围
image(map_tif,col=coltable,xlim=c(75, 130),ylim=c(10,55))
axis(1)#绘制x轴
axis(2)
box()# box() 函数来添加一个方框

plot(map_tif,col=coltable, axes = TRUE,xlim=c(75, 130),ylim=c(10,55))
if(FALSE){
  "
  在 R 语言中，plot() 和 image() 函数都可以用于绘制二维图形，但它们有几个重要的区别：
绘图类型不同：plot() 主要用于绘制散点图、折线图和函数图像等，而 image() 主要用于绘制栅格数据的图像。
数据格式不同：plot() 的输入数据通常为一列或多列数据（如 x 和 y 坐标），而 image() 的输入数据通常为矩阵或数组形式的数据。
绘图方式不同：plot() 默认使用点、线条、文本等来表示数据，而 image() 默认使用颜色来表示数据，并且通常会包含一个色带来显示数据值对应的颜色。
坐标轴设置不同：plot() 可以自由设置坐标轴的范围、标签和刻度等属性，而 image() 的坐标轴设置相对较少，通常只能设置坐标轴标签和刻度。"
}

在这里插入图片描述

#去掉col
plot(map_tif,axes = TRUE,xlim=c(75, 130),ylim=c(10,55))

在这里插入图片描述

#****************************************#
r <- raster("wsiearth.tif")
#去掉col=coltable
plot(r,axes = TRUE,xlim=c(75, 130),ylim=c(10,55))

在这里插入图片描述

library(raster)
map_ras <- raster("wsiearth.tif")
plot(map_ras,col=coltable)
plot(map_ras,col=coltable,xlim=c(75, 130),ylim=c(10,55))
plot(map_ras,col=coltable,xlim=c(75, 130),ylim=c(10,55), maxpixels=5e8)

在这里插入图片描述

# 保存影像
jpeg("../../world.jpg", width=29.7, height=21, units ="cm", res=600) #输入秃瓢
dsn <- paste(wdir,"/data/GlobeData", sep="")
world <- rgdal::readOGR (dsn=dsn, layer="continent")
plot(map_ras, col=coltable, maxpixels=50200200)
plot(world, cex=0.5, add=T)
dev.off()

在这里插入图片描述

1.2.2 GeoTIFF格式

library(rgdal)
library(sp)
setwd(paste(wdir,"/data/wc2", sep=""))
temp01 <- readGDAL("wc2.0_10m_tmin_01.tif")
plot(temp01)

在这里插入图片描述

#col颜色映射，bpy.colors(100) 表示使用 bpy 包提供的颜色映射表，包含 100 种颜色
image(temp01, col=bpy.colors(100))

在这里插入图片描述

#spplot() 函数可以将空间数据图绘制为网格图或条件图即出现条带并非段
spplot(temp01, 
       xlab = "Longitude",
       ylab = "Latitude",  
       par.settings = list(fontsize = list(text = 8)), 
       scales = list(at=seq(-180, 180, 50), draw=TRUE), 
       col.regions=rev(heat.colors(10))[1:9],
       maxpixels=1000000)
#修改色带时，由于函数heat.colors()返回的是从红到黄再到白的颜色渐变，因此我们使用rev()函数翻转这个向量，并且只选择其中的前9个颜色（去掉白色），以得到从蓝到青到绿的颜色渐变

![在这里插入图片描述](https://img-blog.csdnimg.cn/4e27c4f916494c9aa00257c7510ac93a.png

temp01_name<- GDAL.open("wc2.0_10m_tmin_01.tif")
#查看元数据
GDALinfo("wc2.0_10m_tmin_01.tif")
system("gdalinfo wc2.0_10m_tmin_01.tif")

#getRasterTable()函数通常是用于读取栅格数据（如遥感影像）并将其转换为表格格式
temp01_table <- getRasterTable(temp01_name)
head(temp01_table)
#观察数据信息
#str()显示对象（如变量、数据框等）的结构信息，包括对象的类别、属性和组成部分等
str(temp01_table)
summary(temp01_table)
#getRasterData() 函数通常用于读取栅格文件的像素数据，并将其存储在一个 R 语言对象中。
temp01_data <- getRasterData(temp01_name)
#View打开一个新的数据视图窗口查看数据
View(temp01_data)

在这里插入图片描述

1.3 坐标系统

1.3.1 坐标系统

?CRS()
#EPSG 是一个包含 EPSG 投影系统定义的数据框，通过执行 make_EPSG() 函数创建
EPSG <- make_EPSG()
#grep() 函数查找所有包含字符串 "china" 的行，并将它们的索引位置存储在变量 i 中
i <- grep("china", EPSG$note, ignore.case=TRUE)
#从 EPSG 数据框中选取第 1 到第 3 行，其实就是选取前三个包含 "china" 的投影系统定义
EPSG[i[1:3],]
#选取前三个包含 "china" 的投影系统定义的第三列（即 EPSG 代码）
EPSG[i[1:3], 3]

在这里插入图片描述

1.3.2 坐标变换

dsn <- paste(wdir,"/data/GlobeData", sep="")
world <- rgdal::readOGR (dsn=dsn, layer="continent")
#读取坐标参考系统
#proj4string(shapefile)
#st_crs(geojson)
#投影坐标系统
projection(world)
proj4string(world)
#[1] "+proj=longlat +datum=WGS84 +no_defs"
if(FALSE){
  "
  +proj=longlat：使用经度（longitude）和纬度（latitude）作为空间坐标系统，也就是常说的地理坐标系。
+datum=WGS84：指定了该数据集使用的水准面或基准面，即WGS 84大地水准面。
+no_defs：告诉计算机不要使用任何默认的转换参数。
这个CRS描述了一个使用WGS 84大地水准面作为基准的地理坐标系。跟投影无关
  "
}
par(mfrow=c(2,1), mar=c(2,2,2,2)) 
plot(world, col="lightblue")

#CRS()函数用于在R中创建新的CRS对象，带有Robinson投影和WGS84基准的新CRS对象被创建
newcrs <- CRS("+proj=robin +datum=WGS84")
#spTransform()函数用于将空间数据集从一个CRS转换到另一个CRS
#****************************************#
world_New <- spTransform(world, newcrs)
plot(world_New, col="lightblue")

在这里插入图片描述
（1）坐标系理解

GIS中那些坐标系
 概念讲解：大地水准面 | 地球椭球体 | 参考椭球体 | 大地基准面 | 地图投影
 标准地图服务系统

（2）R语言实现
R语言之空间数据操作
 R语言空间数据分析学习笔记5——栅格数据处理

三、空间数据操作

library(rgdal)
library(sf)
library(rgeos)
library(raster)

3.1 矢量数据

3.1.1分割&合并

dsn <- paste(wdir,"/data/china", sep="")
#p <- rgdal::readOGR (dsn=dsn, layer="bou2_4p") #读取.shp文件
#GBK 是一种针对中文的字符编码方式，GBK 编码采用双字节表示一个汉字更占内存
p <- rgdal::readOGR (dsn=dsn, layer="bou2_4p", encoding='GBK') #读取.shp文件
par(mfrow=c(2,2),mar=c(1,1,1,1))
plot(p)
geom(p)

### 分割####
sd <- subset(p, NAME == "山东省")
plot(sd)
js <- subset(p, NAME == "江苏省")
plot(js)

##合并####
sj <- bind(sd,js)
plot(sj, col="red")

在这里插入图片描述

3.1.2 聚合

#按照 NAME 属性进行分组并通过aggrate取并集
#****************************************#
p1 <- aggregate(p, by = "NAME")

plot(p)
plot(p1)
plot(p, col=rainbow(length(p)))
#rainbow(length(p))会生成一个长度(要素个数)为p的向量，其中每个元素都是一个RGB颜色值
plot(p1, col=rainbow(length(p1)))

在这里插入图片描述

3.1.3计算中心点

cen1 <- gCentroid(p1,byid=T)#计算多边形中心点，byid = T 表示按要素 ID 进行分组计算
plot(p1)
points(cen1, pch=16, col="red",cex=0.5)#绘制点图cex是0.5倍大小
##p1@data$NAME矢量对象 p1 中的每个区域的 NAME 属性值作为标注文本
text(cen1,label=p1@data$NAME, adj=0.5, pos=1, cex=0.5)#pos=1在下方，2在左

在这里插入图片描述

3.1.4 切割（1）

#raster 函数将 p 转换成了一个 2 行 2 列、值为 1 到 4 的栅格数据集
par(mfrow=c(2,1),mar=c(2,2,2,2))
z <- raster(p, nrow=2, ncol=2, vals=1:4)
names(z) <- 'Zone'
#as 函数将该栅格数据集转换成了一个空间多边形数据框。每个多边形表示一个区域，并包含与该区域相关的属性信息
z <- as(z, 'SpatialPolygonsDataFrame')
#calss:SpatialPolygonsDataFrame/features:4/names:zone
z
#选择第二个feature即第二个矩形范围
z2 <- z[2,]
plot(p)
plot(z, add=TRUE, border='blue', lwd=5)#z是以范围的矩形划分4块
plot(z2, add=TRUE, border='red', lwd=2, density=3, col='red')
#intersect取交集
p2 <- intersect(p, z2)
plot(p2)#p2是p在第二范围矩形的区域

在这里插入图片描述

3.1.4 切割（2）

选择范围矩形
e <- extent(110, 130, 30, 50)
#****************************************#
#crop裁减 从p中截取e,也可用于栅格
pe <- crop(p, e)
plot(p)
plot(e, col="red", add=T)#在中国地图上画出选择范围显示红框
plot(pe)

在这里插入图片描述

3.2 栅格数据

四、绘制

4.1 plot


library(sp)
library(rgdal)
library(sf)
library(rgeos)
library(raster)

par(mfrow=c(2,2),mar=c(4,2,4,2))
dsn <- paste(wdir,"/data/china", sep="")
p <- rgdal::readOGR (dsn=dsn, layer="bou2_4p") #读取.shp文件
p1 <- aggregate(p, by = "NAME")
#查看字段：[1] "AREA""PERIMETER" "BOU2_4M_" "BOU2_4M_ID" "ADCODE93" "ADCODE99" "NAME"
names(p)
p$PERIMETER
#\xba\xda\xc1\xfa\xbd\xadʡ编码有误修改
#iconv:convert a character vector between encodings: the ‘i’ stands for ‘internationalization’.
P02 <- iconv(p$NAME, from = "GBK", to = "UTF-8")
head(P02)#"黑龙江省" "内蒙古自治区" "新疆维吾尔自治区" "吉林省""辽宁省"           "甘肃省


plot(p1)
plot(p1, col=cm.colors(length(p1)), axes=T)# 有边框和刻度(自调整需要)

plot(p1, col=cm.colors(length(p1)), axes=T, border="blue",)
#(4)
#****************************************#
plot(p1, col=cm.colors(length(p1)), axes=T, border="blue",
     xlim = c(100, 120), ylim = c(25, 45), yaxs="i",xaxs="i",asp=T,main="中国局部",sub="(4)") #作图范围控制在北纬25度至45度,东经100至120度
box()
abline(h=40)#在北纬40度加一条纬线
abline(v=c(105,110,115), col="grey")#abline画线，v指垂直线
abline(h=c(30,35,40), col="grey") 
points(110, 40, col="red", pch= 0) #在东经110度，北纬40度上添加一个红色方块
text(p1, 'NAME', col="black",cex=1)#编码出现问题
help("plot")

在这里插入图片描述

4.2 ggplot

4.2.1 转换为矩阵

library(ggplot2)
library(maptools)
library(plyr)
library(rgdal)
##数据读取####
dsn <- paste(wdir,"/data/china", sep="")
map <- rgdal::readOGR (dsn=dsn, layer="bou2_4p") #读取.shp文件

class(map)#"SpatialPolygonsDataFrame"
##转换####
#SpatialPolygonsDataFrame 类型的对象 map 转换成一个数据框对象，其中每个多边形会被拆分成多个点并按顺序连接起来
names(map)
p <- fortify(map)#空间信息
#
# 查看属性表

```ruby
#**********************#
# head(st_drop_geometry(shapefile)) or
head(map@data)#但是由于name存在编码问题
head(p)
#long      lat order  hole piece id group
#1 121.4884 53.33265     1 FALSE     1  0   0.1

#提取出 map 对象中的属性数据，保存在变量 x 中。
x <- map@data
head(x)
x <- transform(x, NAME = iconv(NAME,from = 'GBK')) #修改name编码问题
#将变量 x 转换为数据框类型，并将原始数据中的行号作为列添加到数据框 xs 中
xs <- data.frame(x,id=rownames(x))
head(xs)
#数据框 p 和数据框 xs 中的行匹配（自动化地根据特定的列名）得到一个新的数据框 my_data，其中包含了空间信息和属性信息
my_data <- join(p, xs, type="full")
head(my_data)
#long      lat order  hole piece id group   AREA PERIMETER BOU2_4M_ BOU2_4M_ID ADCODE93 ADCODE99     NAME
#1 121.4884 53.33265     1 FALSE     1  0   0.1 54.447    68.489      1         23   230000   230000 黑龙江省
#转换结束
class(my_data)#"data.frame"

在这里插入图片描述

4.2.2 投影变换

4.2.3 地图绘制

p <- ggplot(data=<输入数据框>, mapping=aes(x=<变量名>,y=<变量名>,<维度>=<变量名>, <.group=group>/<fill=factor=value.>)) +
#可以省略data =, mapping =, x =, y =
 geom_<图形类型()geom_point/polygon/>(fill="white", colour = "blue") + 
  scale_<映射>_<类型>(<...>) +
  coord_<类型>(<...>) +
  labs(<...>)

##geom_polygon #x = lat, y = long是my_data的信息，group 则指示了哪些点或线应该被视为同一组，以便正确绘制多边形。
ggplot(my_data, aes(y = lat, x = long,group = group))+ geom_polygon(fill="white", colour = "blue")
#二者的区别
#将地图对象转换为数据框，提取其经纬度坐标以及一些额外属性（例如NAME，id等），并将其与空间信息组合。
#这个转换后的数据框存储在my_data对象中。然后，使用这个新的数据框来绘制蓝色的路径，其经纬度坐标用于在对象之间绘制线条。

在这里插入图片描述

4.2.4 保存

最后导出图像即可。ggsvae 提供了丰富的参数定义输出的图像。对于需要投稿 SCI 的文章，通常 Author Guidelines 要求提供不低于 300 DPI 的图片文件。如果允许，保存为 PDF 文件会是不错的方法，毕竟通用性和文件大小都能得到很好的满足。


ggsave(filename = "gmap.pdf", plot = gmap, width = 9 * scale_parm, height = 6.2 * 
    scale_parm, units = "cm", dpi = 600)