利用raster包进行栅格数据处理（2）

,欢迎也在微信上查看相关内容。

数据下载地址：

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/srtm_pnpg.zip

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/BOUNDS_PNPG.zip

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/CIVPARISH_PNPG.zip

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/legend_clc.csv

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/CLC_06_12.zip

library(sp)

library(rgdal)

library(raster)

#数据解压在以下目录中

setwd("F:/test")

rst <- raster("./srtm_pnpg/srtm_pnpg.tif")

maskLayer <- as(readOGR(dsn = "./srtm_pnpg", layer = "pnpg_bounds"), "SpatialPolygons")

plot(rst, main="Elevation (meters) for Peneda-Geres\n National Park", xlab="X-coordinates", ylab = "Y-coordinates")

plot(maskLayer, add=TRUE)

#裁剪

rstMasked <- mask(rst, maskLayer)

plot(rstMasked, main="Elevation (meters) for Peneda-Geres\n National Park",

xlab="X-coordinates",  ylab = "Y-coordinates")

plot(maskLayer, add=TRUE)

栅格聚合是通过将矩形区域中的单元值分组以创建更大/更粗的单元来创建新的栅格层的过程。这种分组可以使用任何用户定义的函数来汇总多值(在矩形区域内)，并提供值(例如，mean, sd, min, max, sum)。Aggregate函数的聚合 “粗度”由fact参数控制，该参数表示每个方向(水平和垂直)的单元数。另外，可以使用两个整数分别表示水平和垂直的聚合因子。

rstAggFact2Mean <- aggregate(rst, fact = 2, fun = mean)

rstAggFact2SD <- aggregate(rst, fact = 2, fun = sd)

rstAggFact7Mean <- aggregate(rst, fact = 7, fun = mean)

rstAggFact7SD <- aggregate(rst, fact = 7, fun = sd)

# Plot the newly aggregated rasters

par(mfrow = c(2,2))

 

plot(rstAggFact2Mean, main = "Aggregation factor = 2 | Mean")

plot(rstAggFact2SD, main = "Aggregation factor = 2 | SD")

 

plot(rstAggFact7Mean, main = "Aggregation factor = 7 | Mean")

plot(rstAggFact7SD, main = "Aggregation factor = 7 | SD")

分区计算

要求输入的数据要具有相同的分辨率和坐标系统。

rstCivPar <- raster("./srtm_pnpg/PNPG_CivilParishes.tif")

par(mfrow = c(1,2))

plot(rst, main="Elevation (meters)")

plot(maskLayer, add=TRUE)

plot(rstCivPar, main="Civil parishes PGNP")

plot(maskLayer, add=TRUE)

交叉表格

执行两个栅格数据集的交叉表非常有用，例如，当想评估两个不同日期之间的土地覆盖变化时。这也是生成混淆矩阵的初步步骤，从混淆矩阵中可以计算出几个分类性能指标。

在这个例子中，我们将使用Corine土地覆盖(CLC)，一个来自欧洲环境局(EEA) 2006年和2012年的数据集，来分析土地覆盖组成的变化。

#交叉表

clcLeg <- read.csv('./srtm_pnpg/legend_clc.csv',stringsAsFactors = F)

# Load the Corine Land cover dataset for 2006 and 2012

clc06 <- raster("./srtm_pnpg/clc2006_100m.tif")

clc12 <- raster("./srtm_pnpg/clc2012_100m.tif")

# 'Ratify' the rasters, i.e., inform that these are

# categorical/discrete datasets

clc06 <- ratify(clc06)

clc12 <- ratify(clc12)

# Perform the crosstab

ct <- crosstab(clc06, clc12, long = TRUE)

# Get the class integer codes and size

lv <- unique(c(ct[,1], ct[,2]))

n <- length(lv)

# Create the square confusion matrix filled with 0's

cm <- matrix(0, nrow = n, ncol = n, dimnames = list(lv,lv))

# Fill the matrix following each line of the contingency table

for(i in 1:nrow(ct)){

  m1<-as.character(ct[i,1])

  m2<-as.character(ct[i,2])

  cm[m1,m2] <- ct[i,3]

}

# Calculate

forestLossAreas <- (clc06 %in% 23:25) & (clc12 %in% c(29,33))

# Plot the results

plot(forestLossAreas, main="Forest loss in PGNP (NW PT)", xlab="x-coord", ylab="y-coord")

plot(spTransform(maskLayer, CRS=crs(clc06)), add=TRUE)