Invest模型获取区域生境质量

Invest模型简介

INVEST模型，全称Integrated Valuation of Ecosystem Services and Tradeoffs（生态系统服务与权衡的综合评估），是由The Natural Capital Project开发的一套用于评估和量化自然资本及其提供的生态系统服务价值的工具。它旨在帮助决策者理解不同土地和海洋使用变化对人类福祉的影响，并为保护、管理和恢复自然资本提供科学依据。

INVEST模型可以支持多种类型的分析，包括但不限于：

• 栖息地质量评估：评估栖息地的质量以及其对物种多样性的支持能力。
• 碳储存与固碳服务：计算特定区域内的碳储量，并预测未来土地利用变化对碳循环的影响。
• 水源涵养功能：评估流域内植被对水文过程的作用，如减少土壤侵蚀、保持水质等。
• 海岸带防护服务：分析沿海湿地和其他自然特征如何保护沿岸地区免受风暴潮和海平面上升的影响。
• 渔业生产潜力：估计海洋或淡水系统中鱼类资源的数量及分布情况。
• 景观美学价值：衡量自然景观对于休闲旅游活动的重要性。

通过这些分析，INVEST能够帮助用户识别关键生态服务供给区，比较不同管理方案的成本效益，从而促进更明智的土地和水资源规划决策。此外，该模型还支持跨学科合作，鼓励科学家、政策制定者、企业和社区成员共同参与自然资源的可持续管理。

Invest模型官方网址：Invest官网
Invest中文文档地址：Invest中文用户指南

Invest模型下载

在Invest模型官方网址中提供了下载链接Invest模型下载链接
Invest模型提供了windows和mac系统的版本。软件下载安装非常简单，这里不再赘述。
我下载的是mac版本的Invest模型，软件打开如下图所示：
可以看到Invest模型可以做很多的分析工作，包括生境质量分析

数据准备与处理

进入生境质量模块，发现需要指定一个工作空间（就是模型运行后输出文件的地方），需要现状土地覆盖的栅格文件（注意：栅格文件的投影方式需要是UTM，是m的投影方式，不是经纬度的投影方式，否则需要进行投影转换）。可选项可以不用准备数据。还需要准备一个威胁表格（同时需要对威胁因子进行像元重分类，即某一威胁因子像元值设为1，其余所有像元值为0）和敏感性表格，格式需要是CSV格式，这两个表格一般需要查找相关文献进行填写。半饱和常数也是查找相关文献进行填写。

了解了需要准备的数据，现在就开始逐一处理。
首先是获取指定区域的土地覆盖数据
我利用GEE平台获取了指定shp区域的土地覆盖数据并导出为tif文件，代码如下：
这里我使用的投影方式是EPSG:4326，是按经纬度投影的，后面会做转换。也可以在GEE中直接实现UTM投影，我所在的区域采用UTM投影应该是采用EPSG:32645

// 加载shp文件
var region = ee.FeatureCollection("projects/ee-pengbr/assets/Kongque10kmbuffer");

// 加载MODIS土地利用数据
var landcover = ee.ImageCollection("MODIS/006/MCD12Q1")
                  .filter(ee.Filter.date('2001-01-01', '2001-12-31'))
                  .first()
                  .select('LC_Type1'); // 选择特定波段;

// 裁剪到shp区域
var landcover_clipped = landcover.clip(region);

// 可视化
var visParams = {
  min: 1,
  max: 17,
  palette: [
    '05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff'
  ]
};

Map.centerObject(region);
Map.addLayer(landcover_clipped, visParams, "Land Cover 2001");

// 导出数据
Export.image.toDrive({
  image: landcover_clipped,
  description: 'LandCover_2001',
  folder: 'GEE',
  region: region.geometry(),
  scale: 250,
  crs: 'EPSG:4326',
  maxPixels: 1e13
});

为了方便后续tif文件的重分类处理及投影转换等，我编写了一个tif文件处理的工具类，上传到了gitee仓库中，地址：tif像元处理工具类
具体代码如下：
包括获取像元平均值、像元总数、像元重分类、shp和tif投影转换（采用python调用linux命令的方式）等方法

import rasterio
import numpy as np
from collections import Counter
from rasterio.features import geometry_mask
import geopandas as gpd
import subprocess


class TifProcess:
    def __init__(self, tif_path, exclude_zero=False):
        """
        初始化类
        :param tif_path: TIFF 文件路径
        :param exclude_zero: 是否排除值为 0 的像元
        """
        self.tif_path = tif_path
        self.exclude_zero = exclude_zero

    # 读入tif文件
    def read_tif(self):
        """
        读取 TIFF 文件
        :return: 数据数组和 nodata 值
        """
        with rasterio.open(self.tif_path) as src:
            data = src.read(1)  # 假设只有一个波段，如果多波段，修改此处
            nodata_value = src.nodata  # 获取无效值标记（如果存在）
            profile = src.profile  # 获取原始文件的元数据
        return data, nodata_value, profile

    # 过滤有效数据
    def _filter_valid_data(self, data, nodata_value):
        """
        过滤掉无效值和零值
        :param data: TIFF 数据数组
        :param nodata_value: 无效值标记
        :return: 过滤后的数据数组
        """
        # 替换无效值为 NaN
        if nodata_value is not None:
            data = np.where(data == nodata_value, np.nan, data)
        # 可选：排除值为 0 的像元
        if self.exclude_zero:
            data = np.where(data == 0, np.nan, data)
        return data

    # 获取有效像元的总数
    def get_valid_pixel_counts(self):
        """
        获取 TIFF 文件有效总像元的个数
        :return: 有效像元数量
        """
        data, nodata_value, _ = self.read_tif()
        data = self._filter_valid_data(data, nodata_value)
        return np.sum(~np.isnan(data))

    # 获取像元平均值
    def get_pixel_mean_value(self):
        """
        获取 TIFF 文件的像元平均值（可排除无效值和 0 值）（针对连续数据，如NDVI值）
        :return: 像元平均值
        """
        data, nodata_value, _ = self.read_tif()
        data = self._filter_valid_data(data, nodata_value)
        return np.nanmean(data)

    # 获取像元种类的个数及占比
    def get_pixel_ratio(self):
        """
        获取 TIFF 文件像元种类的个数、每类的具体个数及占比（针对离散数据，如土地利用类型）
        :return:
            - 像元种类的数量（int）
            - 每类像元的具体个数（dict，键为类别，值为像元数）
            - 每类像元的占比（dict，键为类别，值为占比）
        """
        data, nodata_value, _ = self.read_tif()
        data = self._filter_valid_data(data, nodata_value)
        data = data[~np.isnan(data)]  # 去掉 NaN 值
        counts = Counter(data.ravel())  # 统计每种类别的像元数
        total_pixels = sum(counts.values())  # 总有效像元数量
        ratios = {value: count / total_pixels for value, count in counts.items()}  # 计算占比
        return len(counts), dict(counts), ratios

    # 将指定值变为 1，其他值变为 0（重分类像元）
    def reclassify_pixels(self, target_value, output_path):
        """
        将指定值变为 1，其他值变为 0
        :param target_value: 需要重分类为 1 的目标值
        :param output_path: 保存重分类结果的文件路径
        """
        data, nodata_value, profile = self.read_tif()

        # 重分类操作
        reclassified_data = np.where(data == target_value, 1, 0)

        # 保存结果到新的 TIFF 文件
        profile.update(dtype=rasterio.uint8)  # 更新数据类型为 uint8
        with rasterio.open(output_path, 'w', **profile) as dst:
            dst.write(reclassified_data.astype(rasterio.uint8), 1)

    # 将指定区域外的像元设置为 NoData
    def set_nodata_outside_shp(self, shp_path, output_path):
        """
        将指定矢量区域外的像元设置为 NoData
        :param shp_path: 矢量文件路径（shapefile）
        :param output_path: 保存结果的 TIFF 文件路径
        """
        data, nodata_value, profile = self.read_tif()

        # 如果 nodata 值尚未定义，设置一个默认值
        if nodata_value is None:
            nodata_value = profile.get('dtype', 'float32') == 'int' and -9999 or np.nan

        # 读取矢量文件
        shapes = gpd.read_file(shp_path)

        # 转换矢量到掩膜
        with rasterio.open(self.tif_path) as src:
            mask = geometry_mask(
                [shape.__geo_interface__ for shape in shapes.geometry],
                transform=src.transform,
                invert=True,
                out_shape=(src.height, src.width)
            )

        # 将掩膜外的像元设置为 NoData
        data[~mask] = nodata_value

        # 保存结果到新的 TIFF 文件
        profile.update(nodata=nodata_value)
        with rasterio.open(output_path, 'w', **profile) as dst:
            dst.write(data, 1)

    # 改变tif或shp文件的投影方式
    def reproject(self, extname, target_crs, output_path, shp_path=None):
        global command
        if extname == 'shp':
            command = ['ogr2ogr', '-t_srs', target_crs, output_path, shp_path]
        elif extname == 'tif':
            command = ['gdalwarp', '-t_srs', target_crs, self.tif_path, output_path]
        try:
            subprocess.run(command, check=True, stdout=subprocess.PIPE,
                           stderr=subprocess.PIPE, text=True)
        except subprocess.CalledProcessError as e:
            raise e

现在利用这个工具类查看2001年此区域土地利用类型及其占比

from utools.Tif_pixel_process import TifProcess

tif_path = "LandCover_2001.tif"
tifProcess = TifProcess(tif_path=tif_path, exclude_zero=True)
print(tifProcess.get_pixel_ratio())

结果显示：
一共有6中土地利用类型。我利用的是MODIS数据，土地利用分类表可以参靠我另外一篇博客：D3.js绘制土地利用变化和弦图
12是耕地，10是草地，13是城市用地，16是荒地，11是永久湿地，17是水体。土地利用类型相对简单。其中耕地、城市用地、荒地可以看做是威胁因子，因此需要对其进行重分类。

(6, {12.0: 19430, 10.0: 32738, 13.0: 2795, 16.0: 187834, 11.0: 1086, 17.0: 13}, {12.0: 0.0796651031587234, 10.0: 0.13422934365467248, 13.0: 0.011459802538787024, 16.0: 0.7701397316889166, 11.0: 0.004452717551743366, 17.0: 5.330140715714895e-05})

对威胁因子重分类：
通过以下代码就获取到了耕地、城市和荒地的重分类tif文件。

from utools.Tif_pixel_process import TifProcess

tif_path = "LandCover_2001.tif"
tifProcess = TifProcess(tif_path=tif_path, exclude_zero=False)

tifProcess.reclassify_pixels(target_value=12, output_path='crop_c.tif')
tifProcess.reclassify_pixels(target_value=13, output_path='urban_c.tif')
tifProcess.reclassify_pixels(target_value=16, output_path='barren_c.tif')

转换tif文件的投影方式：

from utools.Tif_pixel_process import TifProcess


def reproject_batching(tif_path, output_path):
    tifProcess = TifProcess(tif_path=tif_path, exclude_zero=False)
    tifProcess.reproject(extname='tif', target_crs='EPSG:32645',
                         output_path=output_path)


reproject_batching('LandCover_2001.tif', 'LandCover_UTM_2001.tif')
reproject_batching('crop_c.tif', 'crops_UTM_c.tif')
reproject_batching('urban_c.tif', 'urban_UTM_c.tif')
reproject_batching('barren_c.tif', 'barren_UTM_c.tif')

接下来就是填写威胁因子表格和敏感性表格了。
威胁因子表如下所示：
在文档中可以查看每个列名所代表的含义，具体数值可以查找相关文献


敏感性表格如下所示：
LULC是土地利用类型编号
habitat (ratio, required): 该LULC类生境的适宜性，0为不适宜，1为完全适宜。
后面三列是威胁因子，[THREAT] (ratio, required): 每个LULC类对每种威胁类型的相对敏感性，其中1表示高敏感性，0表示不受影响。在“威胁表”的“威胁”列中，每个威胁必须有一个对应列。
即使LULC不被认为是生境，也不要将其对每种威胁的敏感性设置为Null或空白，而是输入0。
具体含义可以查看文档

运行模型及结果展示

现在所需要的最基本的数据已经准备好了，接下来就是在Invest软件中填入相关路径运行模型。
半饱和常数这里填入0.5。具体需要查找相关文献和文档。

点击运行后，就会在output文件夹中输出相应文件
其中quality_c.tif表示当前景观的生境质量，deg_sum_c.tif表示当前景观类型的相对生境退化水平。具体请查阅文档。

由于我指定的区域中有大部分生境质量为0的像元，同时在区域外的外接矩阵像元也是显示的0值。为了展示方便，需要将指定shp区域外的值设置为nodata，即透明。可以方便的调用上述的工具类实现：
注意：shp文件的投影方式需要与quality_c.tif投影方式保持一致，这里是UTM投影。经纬度投影方式的shp可以利用上面的工具类进行转换。

from utools.Tif_pixel_process import TifProcess

tif_path = "quality_c.tif"
tifProcess = TifProcess(tif_path=tif_path, exclude_zero=False)
tifProcess.set_nodata_outside_shp(shp_path="孔雀河10km缓冲区_utm.shp",
                                  output_path='quality_2001_c.tif')

最终结果展示：

补充：Invest模型python API

Invest模型提供了python API，通过调用python API实现模型的运行。
接口文档地址：Invest python API

首先安装Invest库，推荐使用conda环境安装。也可以使用pip安装，但是M芯片的Mac电脑不支持。

conda install -c conda-forge natcap.invest

在执行函数中参数是字典格式的。运行脚本与运行软件效果一样。

from natcap.invest.habitat_quality import execute

execute(args={'workspace_dir': './output',
              'lulc_cur_path': 'LandCover_UTM_2001.tif',
              'threats_table_path': 'threats_kongqueriver.csv',
              'sensitivity_table_path': 'sensitivity_kongquerier.csv',
              'half_saturation_constant': 0.5})