Fragstats4.2景观参数计算与景观生态学意义

“ Fragstats4.2是目前最新的生态学参数计算软件，可以很方便的计算景观生态学相关参数，本文以Fragstats4.2介绍生态参数计算及其意义，还有R语言对计算结果进行制图可视化”

01

景观指数计算

—

Fragstats4.2可以和ArcGIS10.0及以下版本进行数据互通，可以直接借助ArcGIS10的Spatial Analysis模块进行Grid数据的读取。然而，ArcGIS10以上版本不再提供Spatial Analysis模块都调用接口，因此10.1及以上的ArcGIS无法再和Fragstats进行ESRI格式的数据交换。尽管GRID格式在高版本ArcGIS时无法和fragstats进行交换通信，Fragstats可以独立使用GeoTIFF格式进行景观格局的分析。
首先新建一个工程，然后Add Layer，如下图所示，导入GeoTIFF格式的数据。

打开Analysis paraments选项卡，在Sampling strategy中选择No sampling, 复选框选择Class metrics、Landscape metrics和Generate patchID file

然后从右侧选择不同的景观参数进行计算。

点击开始计算按钮开始计算，然后再点一下Proceed按钮，这样fragstats就开始计算了

切换到Result面板查看计算结果，然后可以对计算结果进行保存。

02

常用的景观参数介绍

—

斑块类型面积(Total Area, CA/TA)。是度量景观组分和计算其他指标的基础, 其值的大小制约着以某种景观类型斑块为聚集地的物种丰度、数量、食物链及其次生种的繁殖等。

式中，ai为某景观斑块面积，n为斑块数。

斑块个数(Number of Patches, NP)。反映景观的空间格局, 经常被用来描述整个景观的异质性, 其值的大小与景观的破碎度有很好的正相关性, 一般NP 越大, 破碎度越高; NP 越低, 破碎度越低。

NP=n

式中，n为某种景观类型总斑块数。

最大斑块指数（LargestPatch Index, LPI），表征了某一类型的最大斑块占整个景观区域面积都比例。

03

运算结果ArcGIS后处理

—

Fragstats也可以计算得到栅格的结果，Fragstats计算结果会存在空值，在Fragstats里面默认将空值设为-999，为了在ArcGIS里面能够正确的进行显示和分析，需要去除空值。

可以使用栅格计算器将-999设为空值，代码如下：

SetNull("shei.tif" ==  - 999,"shei.tif")

需要注意的是，代码函数和符号应该使用鼠标点击，尤其是-999的负号，如果自己使用键盘输入会报错ERROR 000539

由于生成的TIFF栅格数据很多，可以使用模型构建器批量处理，模型搭建情况如下，有关模型构建器循环处理栅格更多信息请参阅ArcGIS批量转投影

04

景观指数R语言可视化

—

如图，绘制一个斑块类型面积的旋风图，需要将数据转化为如下形式：

TYPE	CA	COUNTRY
其它	157.95	CN
草地	45633.78	CN
其它	160360.56	RU
建设用地	165136.32	CN
建设用地	227573.55	RU
水体	271418.04	CN
湿地	295790.13	CN
草地	575991.81	RU
水体	589802.31	RU
耕地	702245.7	RU
湿地	1016431.74	RU
林地	1832953.86	CN
耕地	3250981.17	CN
林地	10242924.66	RU

使用R语言ggplot2进行绘图：

library(ggplot2)
library(readxl)
CAdata <- read_excel(path = "CA.xlsx", sheet = 2, col_names = T)
p1 <- ggplot(CAdata, aes(x = TYPE,
                   y = ifelse(COUNTRY== 'RU', CA, -CA),
                   fill=COUNTRY))+
  geom_bar(stat = "identity") +
  coord_flip()
  
p1 + scale_y_continuous(limits = c(-11000000, 11000000),
                        breaks = seq(-11000000, 11000000, 10000000),
                        labels = abs)+ 
  labs(y = 'CA(ha)')

p1是未进行坐标轴对称的形式：

完整的代码是修整坐标轴后。