python读取grd数据_Matlab 读取.grd文件

% 【摘要】读取.grd文件；对不同分辨率数据网格插值；对不重合区域取公共交集区

%（1）读 .grd 文件

clc; clear; close all

flistobs = dir('D:\数据\UNH-GRDC\World Runoff Data\runoff_grd\*.grd');

%（2）计算月平均值

for i =2:13 % 留意读取的文件序号

filename = flistobs(i).name;

delimiterIn = ' '; % 设置定界符，从grd文件复制

headerlinesIn = 6; % 设置需要跳过的行数，跳过6行，从第7行开始

datastrct = importdata(filename,delimiterIn,headerlinesIn);

rof = datastrct.data;

rof = flip(rof); % 数据可能是倒置的，可以通过pcolor查看

rof(rof == -9999) = NaN; % 替换缺省值为 Nan

cmprof_mo{i-1} = rof; % 以cell格式存储

end

%（3）计算多年平均值

temp = importdata(flistobs(1).name,delimiterIn,headerlinesIn) % 注意文件在文件夹中的序号

temp.data(temp.data == -9999) = nan;

cmprof_yr = flip(temp.data);

%（4）读经纬度坐标

fpath = 'D:\Models_Output\Clm45Sp\lnd\hist\';

flistsim = dir([fpath,'*.nc']);

fname = flistsim(1).name;

xx = ncread([fpath,fname],'lon');

yy = ncread([fpath,fname],'lat');

for m = 1:length(xx)

if xx(m) >= 180

xx(m) = xx(m) - 360;

else

continue

end

end

xx = [xx(145:end);xx(1:144)]; % 将经度后半段(>=180)拼接到前面

[lonsim,latsim] = meshgrid(xx,yy);

%（5）空间插值（最邻近插值法）

x = -180:0.5:179.5;

y = -55.5:0.5:82.5;

[lonobs,latobs] = meshgrid(x,y);

upcmprof_yrly = griddata(lonobs(:),latobs(:),cmprof_yr(:),lonsim,latsim,'nearest');

%（6）导入模拟值数据

load('simres.mat')

%（7）画图

subplot(2,2,1)

h1 = pcolor(rofsim_yrly)

h1.LineStyle = 'none';

colorbar;

subplot(2,2,2)

h2 = pcolor(upcmprof_yrly)

h2.LineStyle = 'none';

colorbar;

subplot(2,2,3)

h3 = pcolor(rofsim_yrly-upcmprof_yrly)

h3.LineStyle = 'none';

colorbar;

subplot(2,2,4)

hist(rofsim_yrly(:) - upcmprof_yrly(:),500);

%（8）计算两者公共区域（利用Nan 与 非Nan 的计算结果为 Nan 确定交集区域）

comreg = rofsim_yrly - upcmprof_yrly;

rf_simul = rofsim_yrly(~isnan(comreg));

rf_obser = upcmprof_yrly(~isnan(comreg));

%（9）计算评价指标

MAE = mean(abs(rf_obser-rf_simul))

RMSE = sqrt(mean((rf_obser-rf_simul).^2))

noindex = find(rf_obser ~= 0);

delta = abs(rf_obser-rf_simul);

MAPE = mean(delta(noindex)./rf_obser(noindex))

R = corr(rf_obser,rf_simul,'type','Pearson')