腾讯云

文章/答案/技术大牛

发布

社区首页 >专栏 >空间数据处理(一)

空间数据处理(一)

火星娃统计

发布于 2021-03-09 08:03:30

1.8K00

代码可运行

文章被收录于专栏：火星娃统计火星娃统计

运行总次数：0

代码可运行

Spatial 数据简介

Vector data

空间数据的主要形式，类型是点、线和多边形。

点：数据结构为坐标对和附带的值，比如一个地点的温度和它附带的信息比如站点

线：线指的是一系列线段组成的结构，比如河流

多边形：为封闭的折线，起始坐标和终点坐标一致

Raster data

栅格数据通常用于表示空间连续现象，如海拔。栅格将世界划分为大小相同的矩形网格，在遥感数据中称为像素，所有这些网格都有一个或多个值(或缺失值)的变量。栅格单元值通常应该代表它所覆盖区域的平均(或大多数)值或者是中心点的值

与矢量数据相比，栅格数据并不显示存储坐标。通过划分范围来确定，从行数和列数来确定每个单元格的分辨率。

空间数据的简单例子

这里的例子是10个气象站的位置和它们每年的降水量

# 导入和下载需要的包
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# install.packages("rgdal",dependencies = T)
# install.packages("sp")
# install.packages("raster")
library(raster)
library(sp)
# 气象站名字模拟，A-J
name <- LETTERS[1:10]
# 经纬度模拟
longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5,
               -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9,
              36.2, 39, 41.6, 36.9)
# 合并经纬度为矩阵
stations <- cbind(longitude, latitude)

# 模拟降水量数据
# runif函数生成随机正态分布数0-1
set.seed(0)
precip <- round((runif(length(latitude))*10)^3)

# 根据降水量设置点的大小
psize <- 1 + precip/500


# 绘制点图
plot(stations, cex=psize, pch=20, col='red', main='Precipitation')
# 气象站点名字添加
text(stations, name, pos=4)
# 添加图例
breaks <- c(100, 250, 500, 1000)
legend.psize <- 1+breaks/500
legend("topright", legend=breaks, pch=20, pt.cex=legend.psize, col='red', bg='gray')

接下来添加线段和几何图形

lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
# 此对象为几何图形的经纬度
x <- cbind(lon, lat)
# 绘制气象站坐标
# stations 为包含经纬度的矩阵
plot(stations, main='Precipitation')
# 绘制不规则几何图形ploygon函数，可以不用首位点相同
polygon(x, col='blue', border='light blue')
# 将气象站使用线段链接
lines(stations, lwd=3, col='red')
# 绘制几何图形的点的位置
points(x, cex=2, pch=20)
# 绘制气象站的点
points(stations, cex=psize, pch=20, col='red', main='Precipitation')

上述的代码只是简单的实现地图绘制的点、线、几何图形的方式。

Vector 数据

在处理矢量数据的时候，为了方便编写函数，因此定义了很多的类，也就是面向对象，这些类被很多包使用，sp包是处理空间数据的包，虽然sf包也在慢慢完善，但是sp仍然是使用最多的包。

sp包引入了许多类的名称包括：对于矢量数据SpatialPoints，SpatialLines，SpatialPolygons分别代表点、线、几何图形。如果需要包含数据，那么对象为SpatialPointsDataFrame，SpatialLinesDataFrame，SpatialPolygonsDataFrame。接下来从头创建一些空间对象。

SpatialPoints

#建立数据框
longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)
创建spatialpoints对象
library(sp)
pts <- SpatialPoints(lonlat)

# 查看类型
class (pts)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
# 查看具体内容
showDefault(pts)
## An object of class "SpatialPoints"
## Slot "coords":
##       longitude latitude
##  [1,]    -116.7     45.3
##  [2,]    -120.4     42.6
##  [3,]    -116.7     38.9
##  [4,]    -113.5     42.1
##  [5,]    -115.5     35.7
##  [6,]    -120.8     38.9
##  [7,]    -119.5     36.2
##  [8,]    -113.7     39.0
##  [9,]    -113.7     41.6
## [10,]    -110.7     36.9
## 
## Slot "bbox":
##              min    max
## longitude -120.8 -110.7
## latitude    35.7   45.3
## 
## Slot "proj4string":
## CRS arguments: NA

可以看到有一个bbox，这是一个框，显示的是整个数据涵盖的经纬度范围，通过这四个值，经纬度两两匹配，可以确定四个角的经纬度。proj4string是坐标投射的算法，这里没有指定，所以为NA。

指定投射方式为+proj=longlat +datum=WGS84

crdref <- CRS('+proj=longlat +datum=WGS84')
pts <- SpatialPoints(lonlat, proj4string=crdref)

# 使用spatialpoint穿件spatialpointdataframe
precipvalue <- runif(nrow(lonlat), min=0, max=100)
df <- data.frame(ID=1:nrow(lonlat), precip=precipvalue)
# 将降水量数据集与spatialpoints合并
ptsdf <- SpatialPointsDataFrame(pts, data=df)
ptsdf

## class       : SpatialPointsDataFrame 
## features    : 10 
## extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 2
## names       : ID,           precip 
## min values  :  1, 6.17862704675645 
## max values  : 10, 99.1906094830483
# 查看内容
showDefault(ptsdf)
## An object of class "SpatialPointsDataFrame"
## Slot "data":
##    ID    precip
## 1   1  6.178627
## 2   2 20.597457
## 3   3 17.655675
## 4   4 68.702285
## 5   5 38.410372
## 6   6 76.984142
## 7   7 49.769924
## 8   8 71.761851
## 9   9 99.190609
## 10 10 38.003518
## 
## Slot "coords.nrs":
## numeric(0)
## 
## Slot "coords":
##       longitude latitude
##  [1,]    -116.7     45.3
##  [2,]    -120.4     42.6
##  [3,]    -116.7     38.9
##  [4,]    -113.5     42.1
##  [5,]    -115.5     35.7
##  [6,]    -120.8     38.9
##  [7,]    -119.5     36.2
##  [8,]    -113.7     39.0
##  [9,]    -113.7     41.6
## [10,]    -110.7     36.9
## 
## Slot "bbox":
##              min    max
## longitude -120.8 -110.7
## latitude    35.7   45.3
## 
## Slot "proj4string":
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

SpatialLines 和 SpatialPolygons

制作线，这里的函数均来自raster

lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(lon, lat)
# 使用splines函数制作线
lns <- spLines(lonlat, crs=crdref)
lns
## class       : SpatialLines 
## features    : 1 
## extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

制作几何图形,几何图形的对象结构较为复杂，不做展开

# 使用函数sppolygons
pols <- spPolygons(lonlat, crs=crdref)
pols
## class       : SpatialPolygons 
## features    : 1 
## extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs
绘制线和几何图形
plot(pols, axes=TRUE, las=1)
plot(pols, border='blue', col='yellow', lwd=3, add=TRUE)
points(pts, col='red', pch=20, cex=3)

Raster 数据

栅格数据处理主要使用的是raster包。raster包主要的三个对象，RasterLayer，RasterBrick，RasterStack。

RasterLayer

RasterLayer对象表示单层栅格数据。一个RasterLayer对象存储一些描述它的基本参数。这些参数包括列和行数、空间范围和坐标参考系统。此外，RasterLayer可以存储单元值的文件的信息。

创建RasterLayer

# 创建一个10行10列的栅格数据框架
r <- raster(ncol=10, nrow=10, xmx=-80, xmn=-150, ymn=20, ymx=60)
r
## class      : RasterLayer 
## dimensions : 10, 10, 100  (nrow, ncol, ncell)
## resolution : 7, 4  (x, y)
## extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs
# 给栅格数据单元添加值
values(r) <- 1:ncell(r)
r
## class      : RasterLayer 
## dimensions : 10, 10, 100  (nrow, ncol, ncell)
## resolution : 7, 4  (x, y)
## extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 100  (min, max)

绘制图形

plot(r)
# 添加几何和线段
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(lon, lat)
# 绘制几何图形
pols <- spPolygons(lonlat, crs='+proj=longlat +datum=WGS84')
# 添加点
points(lonlat, col='red', pch=20, cex=3)

# 绘制几何点
plot(pols, border='blue', lwd=2, add=TRUE)

RasterStack和RasterBrick

在大多数的情况下，使用的是单层的栅格数据分析，但是在一些案例中，需要使用到多层数据，因此引入RasterStack和RasterBrick。RasterStack针对的是单一的多层文件，RasterBrick针对的是多个文件

事实上，Rasterstack是具有相同空间范围和分辨率的RasterLayer对象的集合。

RasterBrick是一个真正的多层对象，处理RasterBrick比处理表示相同数据的Rasterstack更有效.

制作RasterStack

# r是rasterlayer
r2 <- r * r
r3  <- sqrt(r)

# 使用stack函数，建立rasterstack对象
s <- stack(r, r2, r3)
s

## class      : RasterStack 
## dimensions : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
## resolution : 7, 4  (x, y)
## extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : layer.1, layer.2, layer.3 
## min values :       1,       1,       1 
## max values :     100,   10000,      10
# 绘制图像
plot(s)

可以看到对象s有三个层,使用brick可以转换stick对象为brick

b <- brick(s)
b
## class      : RasterBrick 
## dimensions : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
## resolution : 7, 4  (x, y)
## extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer.1, layer.2, layer.3 
## min values :       1,       1,       1 
## max values :     100,   10000,      10