InVEST 模型相关实验

我真的真的真的会喜欢读博

已于 2024-10-02 22:09:46 修改

阅读量452

点赞数 7

分类专栏：数据处理生态参数反演文章标签： arcgis

于 2024-09-21 16:43:05 首次发布

本文链接：https://blog.csdn.net/RS_LAB/article/details/142418912

版权

数据处理同时被 2 个专栏收录

3 篇文章 3 订阅

订阅专栏

生态参数反演

1 篇文章 1 订阅

订阅专栏

博文仅记录、分享自己的实验步骤，如有错误还请不吝赐教。

1. 产水量模块

1.01 工作空间

1.02 批处理后缀文件

1.03 降水数据

（1）1901-2023年中国1km分辨率逐月降水量数据集，源于国家地球系统科学数据中心

http://www.geodata.cn/main/face_science_detail?guid=192891852410344&typeName=face_science

原始数据为CRU，降水量单位为0.1mm。

（2）中国区域1km分辨率逐月近地表平均气温数据集（2000-2020年），源于国家地球系统科学数据中心

原始数据为ERA5

https://www.geodata.cn/main/face_science_detail?id=56244&guid=3034046

1.04 蒸散发数据

(1) 1901-2023年中国1km分辨率逐月潜在蒸散发数据集，源于国家地球系统科学数据中心

空间约1km，时间为1901.1-2023.12，单位为0.1mm。

http://www.geodata.cn/main/face_science_detail?guid=34595274939620&typeName=face_science

1.03 降水数据、1.04 蒸散发数据部分数据处理：

参考：（1）ArcGIS打开一个nc文件，使用多维工具随便提取一个样例数据

将单波段的结果输出为tiff影像

下面的代码需要修改的有：

工作空间：workingPath，存放nc文件的文件夹。

样例数据：sample，你刚刚提取出的单波段tif影像路径

变量：例如气温改为‘tmp’，降水改为‘pre’，潜在蒸散发改为‘etp’等等。

比例系数：例如0.1mm，则所有像元×0.1，（etp=etp*0.1;）

（可以从GIS中多维工具查看nc文件的变量是什么）

clc;
clear;
workingPath = 'E:\PET\AOrigion';% 设置工作路径

sample='E:\PET\Sample2.tif'; %%刚刚提取的样例数据的地址，用于读取其坐标系统
[A,R]=geotiffread(sample);
info=geotiffinfo(sample);
[m,n]=size(A);
cd(workingPath);
files = dir('pet_*.nc');

% 创建输出文件夹 'output'
outputFolder = fullfile(workingPath, 'output');
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end


for k = 1:length(files)  % 循环读取每个.nc文件
    [~, name, ~] = fileparts(files(k).name);
    year = str2double(name(5:8));    % 获取文件名和年份参考 "pet_2014.nc"
  
    % 读取.nc文件中的数据
    ncfile = fullfile(workingPath, files(k).name);
    lon = ncread(ncfile, 'lon');
    lat = ncread(ncfile, 'lat');
    time = ncread(ncfile, 'time');
    etp = ncread(ncfile, 'etp');  %这里需要修改降水修改为pre，温度为tmp，GIS多维工具中查看
    etp(etp == -32768) = NaN;
    etp= etp*0.1;%注意比例系数，如果没有则可以不填

    datasum=zeros(m*n,length(time))+NaN; 

    for m = 1:12
        % 提取每个月的数据
        if m <= length(time)
            monthly_data = etp(:, :, m)';            
            % 设置输出文件名
            filename = sprintf('etp_%d_month_%02d.tif', year, m);
            outputFile = fullfile(outputFolder, filename);

            geotiffwrite(outputFile,monthly_data,R,'TiffType','bigtiff','GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)      
        end
    end
end

1.05 Root restricting layer depth（根系限制层）

数据所需要的单位为毫米（mm），注意使用数据的单位转化（栅格计算器）。

（1）SoilGrids250m 2017-03 - Depth to bedrock (R horizon) 基岩深度

https://data.isric.org/geonetwork/srv/eng/catalog.search#/metadata/bfb01655-db81-4571-b6eb-3caae86c037a

（2）Depth-to-bedrock map of China at a spatial resolution of 100 meters （中国基岩深度，100m空间分辨率，数据单位为m，需要×1000）

Depth-to-bedrock map of China at a spatial resolution of 100 meters | Scientific DataDepth to bedrock influences or controls many of the Earth’s physical and chemical processes. It plays important roles in soil science, geology, hydrology, land surface processes, civil engineering, and other related fields. However, information about depth to bedrock in China is very deficient, and there is no independent map of depth to bedrock in China currently. This paper describes the materials and methods to produce high-resolution (100 m) depth-to-bedrock maps of China. For different research and application needs, two sets of data are provided for users. One is the prediction by the ensemble of the random forests and gradient boosting tree models, and the other is the prediction and the uncertainty of prediction based on quantile regression forests model. In comparison with depth-to-bedrock maps of China extracted from previous global predictions, our predictions showed higher accuracy and more spatial details. These data sets can provide more accurate information for Earth system research compared with previous depth-to-bedrock maps. Machine-accessible metadata file describing the reported data: https://doi.org/10.6084/m9.figshare.11358929https://www.nature.com/articles/s41597-019-0345-6

Yan, Fapeng; Shangguan, Wei; Zhang,, Jing; Hu, Bifeng (2019). Depth-to-bedrock map of China at a spatial resolution of 100 meters. figshare. Collection. https://doi.org/10.6084/m9.figshare.c.4714514.v1

对应的论文：

Yan, Fapeng, et al. "Depth-to-bedrock map of China at a spatial resolution of 100 meters." Scientific Data 7.1 (2020): 2.

1.06 Plant Available Water Content（植物可利用含水量）

（1）官方文档计算

找到 SoilGrids250m 2017-03 - Derived available soil water capacity (volumetric fraction) until wilting point，

1）下载（无论怎么样下载都很慢，我下了整整一天）

2）裁剪：根据使用手册，将其重新命名，按规则命名。一定要建立一定区域的缓冲区进行掩膜提取。

裁剪输出的名称为： AWC_sl1_clip.tif, AWC_sl2_clip.tif … AWC_sl7_clip.tif.

3 ）计算

将裁剪后的文件全部加入图层。打开栅格计算器，输入公式，输出文件名称等等。

桌面版ArcGIS公式（InVEST使用手册）：

0.005 * 0.5 * (5 * ("AWC_sl1_clip.tif"+"AWC_sl2_clip.tif")+10 * ("AWC_sl2_clip.tif"+"AWC_sl3_clip.tif")+15*("AWC_sl3_clip.tif"+"AWC_sl4_clip.tif")+30 * ("AWC_sl4_clip.tif" + "AWC_sl5_clip.tif")+40 * ("AWC_sl5_clip.tif"+"AWC_sl6_clip.tif") + 100 * ("AWC_sl6_clip.tif" + "AWC_sl7_clip.tif"))

4 ）使其适配InVEST模型

计算出来的结果应该在0~100之间（检查），并使用栅格计算器转换为百分比

后续的操作使其适配土地利用数据包括以下步骤：

重投影（例如UTM WGS84 48N 坐标系） ——重采样（例如双线性采样至30m）————掩膜提取（注意设置处理范围，使其和土地利用数据行列号一致）

（2）经验公式估算

周文佐.基于GIS的我国主要土壤类型土壤有效含水量研究[D].南京农业大学,2003.

周文佐,刘高焕,潘剑君.土壤有效含水量的经验估算研究——以东北黑土为例[J].干旱区资源与环境,2003,(04):88-95.

Zhou, Wenzuo, et al. "Distribution of available soil water capacity in China." Journal of Geographical Sciences 15 (2005): 3-12.

Available soil water capacity (ASWC) is the volume of soil water that should be available to
plants. （Zhou, Wenzuo, et al. "Distribution of available soil water capacity in China." Journal of Geographical Sciences 15 (2005): 3-12.）定义差别不大。

图片摘自（Zhou, Wenzuo, et al. "Distribution of available soil water capacity in China." Journal of Geographical Sciences 15 (2005): 3-12.）中国ASWC处于30以下。

图片摘自（Zhou, Wenzuo, et al. "Distribution of available soil water capacity in China." Journal of Geographical Sciences 15 (2005): 3-12.）。有两个公式可以计算ASWC。

所需数据为土壤数据库中的CLAY、SAND、SILT、OC。

由于我的研究区内S层土壤缺少有机碳数据，所以我只取表层土壤计算（HWSD 1.2）。（如果两层数据均有（或采用了HWSD 2.0的7层数据），建议加权求和，例如表层*0.5+深层*0.5）


WORKPath = 'F:\Gannan\Gannan_auxiliary\HWSD1_Gannan\ccopy\';
outputfilename=[WORKPath,'PAWC_T.TIF'];

A_SAND=[WORKPath,'HWSD1_T_SAND.tif'];
B_SLIT=[WORKPath,'HWSD1_T_SILT.tif'];
C_CLAY=[WORKPath,'HWSD1_T_CLAY.tif'];
D_OM=[WORKPath,'HWSD1_T_OM_oc1.724.tif'];
% D_OC=[WORKPath,'HWSD1_T_OC.tif'];%如果你是OC有机碳数据，有机质含量0M=OC*1.724
% [OC, R] = geotiffread(D_OC);
% OM=OC*1.724;

nandata=-128;
[SAND, R] = geotiffread(A_SAND);
SAND(SAND==nandata)=0;
[SLIT, R] = geotiffread(B_SLIT);
SLIT(SLIT==nandata)=0;
[CLAY, R] = geotiffread(C_CLAY); 
CLAY(CLAY==nandata)=0;
[OM, R] = geotiffread(D_OM);
OM(OM<0)=0;

PAWC_CAL=54.509-0.132.*SAND-0.003.*SAND.*SAND-0.055*SLIT-0.006.*SLIT.*SLIT ...
-0.738.*CLAY+0.007.*CLAY.*CLAY-2.688.*OM+0.501.*OM.*OM;

geotiffwrite(outputfilename,PAWC_CAL,R)

PAWC计算结果如下：部分区域值特别高，因为有机质含量很高，此时可以利用栅格计算器的Con（）函数进行处理，大于30的部分直接赋予30，或者×0.0001。

经验公式计算出的正常值与InVEST官方文档公式计算结果差别很小，但是经验公式计算出的异常值需要进行修正（因为值的范围在0-100之间）。

注：数据可能存在空值

制作一个全是0的图层（研究区shp转栅格，重分类）

找到数据管理工具中的镶嵌至新栅格

镶嵌属性中选择Fist，此时一定要是PAWC在上方：

填充完毕：

1.07 土地利用数据不讲了

1.08 生物物理量表

（1）表的样例

查看手册，产水量模块需要以下几种数据

打开表格

绿色区域为所需要的数据

1.09 Z参数

经验参数，需要根据实际情况慢慢调整（看别人的论文）。

1.10 流域数据

可以使用研究区边界，将县界进行融合，或者使用市界，并在其属性表中加入“ws_id”整型字段。

1.11 子流域数据

可以使用SWAT模型进行子流域划分。

1.2 RUN

检查数据：

年总降水

年总潜在蒸散发

根系限制层

植物可利用含水量

土地利用

参数表、Z参数、

1.3 Validation

验证一般使用水资源公报，可以制作这样一个excel表，对Z参数进行测试

1.4 水源涵养模块

1.4.1 土壤饱和导水率

利用SPAW软件，下载地址：

https://www.nrcs.usda.gov/resources/tech-tools/spaw-version-602

数据处理：利用掩膜提取

参考文章：CSDNhttps://mp.csdn.net/mp_blog/creation/editor/142418912

将BIL数据进行裁剪，生成小区域的tif文件

生成唯一值后输出，后期要进行查找表

将这两个表相关联，可以用vlookup，我是用MATLAB

得到了一个新表，将SAND和CLAY输入SPAWT软件中，逐一计算土壤饱和导水率

查找表获得土壤饱和导水率

1.5 地形湿度指数（topographic wetness index，TWI）

（1）填洼

输入DEM，参数默认

（2）流向

输入填洼后的DEM，参数默认

（3）流量

输入流向，参数默认

（4）栅格计算器

TWI=ln[SCA / tan(slope)]

计算“单位等高线长度集水面积”（SCA，单元栅格汇流面积）和“百分比坡度”。

2.0 水土流失模块

2.1 工作空间

2.2 批处理后缀文件名

2.3 Digital Elevation Model（数字高程模型）

2.4 Erosivity（降雨侵蚀力指数）

章文波,付金生.不同类型雨量资料估算降雨侵蚀力[J].资源科学,2003,(01):35-41.

2.5 Soil Erodibility（土壤侵蚀力指数）

（图片源于：Williams, J. R., and J. G. Arnold. "A system of erosion—sediment yield models." Soil technology 11.1 (1997): 43-55. ）

ArcGIS桌面版公式（请注意栅格文件名）：

先算Soil_Erodibility_K_epic.tif：

(0.2+0.3*Exp(-0.0256*"HWSD1_T_SAND.tif"*(1-"HWSD1_T_SILT.tif"/100)))*Power(("HWSD1_T_SILT.tif"/("HWSD1_T_CLAY.tif"+"HWSD1_T_SILT.tif")),0.3)*(1.0-((0.25*"HWSD1_T_OC.tif")/("HWSD1_T_OC.tif"+Exp(3.72-2.95*"HWSD1_T_OC.tif"))))*(1.0 - ((0.7*(1-"HWSD1_T_SAND.tif"/100))/((1-"HWSD1_T_SAND.tif"/100)+Exp(22.9*(1-"HWSD1_T_SAND.tif"/100)-5.51))))

再算K：

(-0.01383+0.51575*"Soil_Erodibility_K_epic.tif")*0.1317

以及MATLAB代码，两种方法我都试过，计算结果差不多一样。

clc;
clear;

WORKPath = 'F:\Gannan\Gannan_auxiliary\HWSD1_Gannan\TRS\';
outputfilename=[WORKPath,'Soil_Erodibility_K.TIF'];

A_SAND=[WORKPath,'HWSD1_T_SAND.tif'];
B_SLIT=[WORKPath,'HWSD1_T_SILT.tif'];
C_CLAY=[WORKPath,'HWSD1_T_CLAY.tif'];
D_OC=[WORKPath,'HWSD1_T_OC.tif'];
%% SAN、SIL、CLA、C 分别表示土壤砂粒百分含量（%）、
%% 粉粒百分含量（%）、粘粒百分含量（%）和表层（0–30cm）土壤有机碳百分含量（%）
%%  SNI 由1–SAN/100 的关系式计算。
%%
nandata=-128;
[SAND, R] = geotiffread(A_SAND);
SAND(SAND==nandata)=0;
[SLIT, R] = geotiffread(B_SLIT);
SLIT(SLIT==nandata)=0;
[CLAY, R] = geotiffread(C_CLAY); 
CLAY(CLAY==nandata)=0;
[OC, R] = geotiffread(D_OC);
OC(OC<0)=0;

SAND=double(SAND);
SLIT=double(SLIT);
CLAY=double(CLAY);
OC=double(OC);
SNI=1-SAND/100;
% 只处理 SAND > 0 的像元
valid_mask = SAND > 0;
% 初始化 K_EPIC 为 NaN，大小与 SAND 一致K_EPIC = NaN(size(SAND));
% K_EPIC 的计算
K_EPIC(valid_mask) = (0.2 + 0.3 .* exp(-0.0256 .* SAND(valid_mask) .* (1 - SLIT(valid_mask) / 100))) ...
    .* (SLIT(valid_mask) ./ (CLAY(valid_mask) + SLIT(valid_mask))).^0.3 ...
    .* (1.0 - (0.25 .* OC(valid_mask)) ./ (OC(valid_mask) + exp(3.72 - 2.95 .* OC(valid_mask)))) ...
    .* (1.0 - (0.7 .* SNI(valid_mask)) ./ (SNI(valid_mask) + exp(22.9 .* SNI(valid_mask) - 5.51)));

K = (-0.01383 + 0.51575 .* K_EPIC) .* 0.1317;
%geotiffwrite(outputfilename, K, R, 'TiffType', 'bigtiff');
geotiffwrite(outputfilename, K, R, 'CoordRefSysCode', 32648);%% https://epsg.io/