发一个Matlab计算PDSI的函数

kingjia · 发表于 2016-4-7 20:19:43

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 kingjia 于 2016-4-9 18:48 编辑

PDSI计算比较专业权威的计算程序一般认为是美国内布拉斯加大学－林肯分校的专业研究人员发布的计算程序，但没有matlab版本，windows版本运行的时候每次应该只能单站点运算。自己写了一个用于多站点循环调用的matlab函数。
该函数主要为了自己用，看有人求用matlab计算的代码，个人测试可以，但限于本人能力，可能还存在问题与不足，发出来给大家提供一个参考。

说明：
1.该程序是自己写的一个便于多站点循环调用的PDSI计算函数，函数并非PDSI底层计算源码，只是在函数内部调用美国内布拉斯加大学－林肯分校的专业研究人员发布的scpdsi.exe程序进行相关计算。
2.函数输入参数可选：4个或6个或7个。
4个输入参数（土壤有效水持水量，纬度，月平均气温，月总降水）时对应只可以设置一个输出，即月尺度的各种指数（sc_PDSI, sc_PHDI, sc_WPLM, sc_ZIND, PDSI, PHDI, WPLM and Z-index）。
6个输入参数（土壤有效水持水量，纬度，月平均气温，月总降水，weekly平均气温，weekly总降水）时对应只可以设置1-3个输出：月尺度的各种指数，1星期尺度的各种指数，和1星期尺度的CMI指数。
7个输入参数（土壤有效水持水量，纬度，月平均气温，月总降水，weekly平均气温，weekly总降水，1或2或4或13星期尺度），可输出月尺度的各种指数，n星期尺度的各种指数，和1星期尺度的CMI指数。
3.降水和气温可选择两种输入格式，函数内自动识别。
格式1：第一列年份，后边2:13（或2:53）列分别是1-12月份（或1-52星期）对应的值，该格式是exe程序默认需要的格式。
格式2：第一列年份，第二列月份（或星期），第三列对应的值。个人感觉用自己的观测数据处理的时候该格式比较容易组织，因此在写的PDSI函数内写了自动转换代码，考虑了缺值数据，即输入矩阵中缺少某年份某月份的值时，函数内部自动补空值。
4.降水单位mm，气温单位℃。exe程序接受inch 和℉，一般用不到，所以函数暂时没有考虑。
5.原exe程序基于C++写的，测试中发现计算CMI时会出现'-1.#J'，'-1.#IND'等值，Matlab直接dlmread读取会出现问题，个人能力有限，暂时加了一个单独的cmiread函数进行读取和转换（内部用了textscan函数先读取数据到cell中，再在cell中替换字符为空，最后转换为矩阵）。在要求输出CMI的时候运算效率会有所下降。

附件是函数调用例子，和两种格式的输入数据，最后是相关的函数和程序，pdsi.m为主函数。

example.m (369 Bytes, 下载次数: 180)

玺驴驴 · 发表于 2018-6-27 10:16:27

http://westdc.westgis.ac.cn/data/611f7d50-b419-4d14-b4dd-4a944b141175
找到了AWC数据的下载网址请大家自行注册下载

kongfeng0824 · 发表于 2016-4-8 21:15:46

谢谢分享

kongfeng0824 · 发表于 2016-4-8 21:15:48

谢谢分享

liyupeng · 发表于 2016-4-9 13:14:49

楼主你好单独求月指标的时候输入如下总是会出错
MonthlyIndices=pdsi(AWC,lati,mt,mp);
出错显示为
Error using rmdir
G:\aaa\out_tmp\weekly is not a directory.

Error in pdsi (line 233)
rmdir('out_tmp\weekly\','s');

kingjia · 发表于 2016-4-9 18:41:05

liyupeng 发表于 2016-4-9 13:14
楼主你好单独求月指标的时候输入如下总是会出错
MonthlyIndices=pdsi(AWC,lati,mt,mp);
出错显示为

把pdsi.m文件中233行：rmdir('out_tmp\weekly\','s');
替换为：
if nargout>=6
rmdir('out_tmp\weekly\','s');
end

liyupeng · 发表于 2016-4-10 00:35:41

kingjia 发表于 2016-4-9 18:41
把pdsi.m文件中233行：rmdir('out_tmp\weekly\','s');
替换为：
if nargout>=6

谢谢楼主的指导我在试一试

liyupeng · 发表于 2016-4-10 00:35:45

kingjia 发表于 2016-4-9 18:41
把pdsi.m文件中233行：rmdir('out_tmp\weekly\','s');
替换为：
if nargout>=6

谢谢楼主的指导我在试一试

kingjia · 发表于 2016-4-10 18:14:35

liyupeng 发表于 2016-4-10 00:35
谢谢楼主的指导我在试一试

有问题尽管反馈，大家一起不断完善

zhchlue · 发表于 2016-4-14 12:47:03

除了楼主常用的部分，其他部分应该存在不少问题，需要进一步检验

zhchlue · 发表于 2016-4-14 12:47:21

当然也希望大家多多提出讨论

		自动登录	找回密码
密码			立即注册

[源程序] 发一个Matlab计算PDSI的函数

登录后查看更多精彩内容~

评分