1 Star 2 Fork 0

damowangyang/Matlab scripts for GLEAM

加入 Gitee
与超过 1200万 开发者一起发现、参与优秀开源项目,私有仓库也完全免费 :)
免费加入
文件
该仓库未声明开源许可证文件(LICENSE),使用请关注具体项目描述及其代码上游依赖。
克隆/下载
ET_Ratio_Calculation_V35.m 6.71 KB
一键复制 编辑 原始数据 按行查看 历史
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% This scripts is for estimate the ratio of GLEAM_Version3.5a/CLM_Output.
% 2021-02-26 by Dayon
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% clear everything
clc
clear all
yearStart = 1986;
yearEnd = 1999;
nYear = yearEnd - yearStart + 1;
monthStart = 1;
monthEnd = 12;
nMonths = monthEnd - monthStart + 1;
if monthStart < 10
monthStartStr = ['0', num2str(monthStart)];
else
monthStartStr = num2str(monthStart);
end
if monthEnd < 10
monthEndStr = ['0', num2str(monthEnd)];
else
monthEndStr = num2str(monthEnd);
end
startEndDate = ['yr_', num2str(yearStart),'_',num2str(yearEnd),'_mo_',monthStartStr,'_',monthEndStr];
ET_ratio_min = 0.01;
ET_ratio_max = 100;
res = 0.25;
% longtude and latitude
lon = 0.5:1:359.5;
lat = -89.5:1:89.5;
nLon = length(lon);
nLat = length(lat);
% ==================================================================================
% extract CLM5.0 Output Data
% ==================================================================================
inDir = '/home/dayon/Downloads/GLEAM/ModelOutput';
temp2 = zeros(nLon, nLat, nMonths);
monthDays = [31 28 31 30 31 30 31 31 30 31 30 31];
for iMonth = monthStart : monthEnd
if iMonth < 10
iMonthStr = ['0', num2str(iMonth)];
else
iMonthStr = num2str(iMonth);
end
ET = zeros(nLon, nLat, nYear);
for iYear = yearStart : yearEnd
iYearStr = num2str(iYear);
iYearInd = iYear - yearStart + 1;
inFile = ['CLM50SpGs_res1D_CTRL.clm2.h0.',iYearStr,'-',iMonthStr,'.nc'];
inDirFile = [inDir, '/', inFile];
temp1 = ncread(inDirFile, 'QVEGE') + ncread(inDirFile, 'QVEGT') ...
+ ncread(inDirFile, 'QSOIL');
temp1(temp1 > 10000) = nan;
% convert units from mm/s to mm/month
ET(:, :, iYearInd) = monthDays(iMonth) * 3600 * 24 * temp1;
end
temp2(:,:,iMonth) = reshape(mean(ET, 3), nLon, nLat);
end
ET_Model_MonthMean = temp2;
clear temp1 temp2
ET_Model_MonthMean_max = max(max(max(ET_Model_MonthMean)));
ET_Model_MonthMean_min = min(min(min(ET_Model_MonthMean)));
% ==================================================================================
% extract GLEAM Observation
% ==================================================================================
inDir = ['/home/dayon/Downloads/GLEAM/res1D_GLEAM'];
inFile = ['E_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
inDirFile = [inDir, '/' , inFile];
% choose the E_Obs Type
ET_Option = 2;
switch ET_Option
case 1
disp('You choose the Actual evaporation (E) as the Observation ET');
ET_Obs = ncread(inDirFile, 'E');
case 2
inFileEt = ['Et_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
inDirFileEt = [inDir, '/' , inFileEt];
inFileEi = ['Ei_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
inDirFileEi = [inDir, '/' , inFileEi];
inFileEb = ['Eb_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
inDirFileEb = [inDir, '/' , inFileEb];
disp('You choose the sum of Transpiration, Interception loss and Bare-soil evaporation (Et + Ei + Eb) as the Observation ET');
ET_Obs = ncread(inDirFileEt, 'Et') + ncread(inDirFileEi, 'Ei') + ncread(inDirFileEb, 'Eb');
end
ET_Obs_MonthMean = zeros(nLon, nLat, nMonths);
for iMonth = monthStart : monthEnd
ipMonStart = (yearStart - 1980)*12 + iMonth;
ipMonEnd = (yearEnd - 1980)*12 + iMonth;
ipArray = [ipMonStart : 12 : ipMonEnd];
temp3 = ET_Obs(:, :, ipArray);
ET_Obs_MonthMean(:, :, iMonth) = mean(temp3, 3);
end
ET_Obs_MonthMean_min = min(min(min(ET_Obs_MonthMean)));
ET_Obs_MonthMean_max = max(max(max(ET_Obs_MonthMean)));
% ==================================================================================
% Ratio
% ==================================================================================
ET_ratio = ET_Obs_MonthMean./ET_Model_MonthMean;
ET_ratio(isnan(ET_ratio)) = 1.0;
ET_ratio(ET_ratio <= 0) = 1.0;
ET_ratio(ET_ratio >= ET_ratio_max) = ET_ratio_max;
ET_ratio(ET_ratio <= ET_ratio_min) = ET_ratio_min;
ET_ratio_max = max(max(max(ET_ratio)));
ET_ratio_min = min(min(min(ET_ratio)));
% ==================================================================================
% Make netcdf
% ==================================================================================
outDir = '/home/dayon/Downloads/GLEAM/ET_ratio';
switch ET_Option
case 1
outFile = ['Global_ETactual_ratio_GLEAM_', startEndDate, '.nc'];
case 2
outFile = ['Global_ETthree_ratio_GLEAM_', startEndDate, '.nc'];
end
outDirFile = [outDir, '/', outFile];
outid = netcdf.create(outDirFile, 'CLOBBER');
% define dimension
latdimID = netcdf.defDim(outid, 'lat', nLat);
londimID = netcdf.defDim(outid, 'lon', nLon);
timedimID = netcdf.defDim(outid, 'time', nMonths);
% define variable
latid = netcdf.defVar(outid, 'lat', 'double', latdimID);
lonid = netcdf.defVar(outid, 'lon', 'double', londimID);
timeid = netcdf.defVar(outid, 'time', 'double', timedimID);
ETid = netcdf.defVar(outid, 'ET_ratio', 'double', [londimID, latdimID, timedimID]);
QVEGEid = netcdf.defVar(outid, 'QVEGE_ratio', 'double', [londimID, latdimID, timedimID]);
QVEGTid = netcdf.defVar(outid, 'QVEGT_ratio', 'double', [londimID, latdimID, timedimID]);
QSOILid = netcdf.defVar(outid, 'QSOIL_ratio', 'double', [londimID, latdimID, timedimID]);
netcdf.endDef(outid);
netcdf.putVar(outid, latid, lat);
netcdf.putVar(outid, lonid, lon);
timeVar = 0:1:11;
netcdf.putVar(outid, timeid, timeVar);
netcdf.putVar(outid, ETid, ET_ratio);
netcdf.putVar(outid, QVEGEid, ET_ratio);
netcdf.putVar(outid, QVEGTid, ET_ratio);
netcdf.putVar(outid, QSOILid, ET_ratio);
netcdf.reDef(outid);
netcdf.putAtt(outid,latid,'long_name','latitude');
netcdf.putAtt(outid,latid,'units','degrees_north');
netcdf.putAtt(outid,latid,'standard_name','latitude');
netcdf.putAtt(outid,lonid,'long_name','longitude');
netcdf.putAtt(outid,lonid,'units','degrees_east');
netcdf.putAtt(outid,lonid,'standard_name','longitude');
netcdf.putAtt(outid,timeid,'long_name','month');
netcdf.putAtt(outid,timeid,'units','months since 1980-01-01 00:00:00');
netcdf.putAtt(outid,timeid,'standard_name','month');
undef = -9999.9;
netcdf.putAtt(outid,ETid,'long_name','ET_ratio');
netcdf.putAtt(outid,ETid,'units','unitless');
netcdf.putAtt(outid,ETid,'_FillValue',undef);
netcdf.putAtt(outid,QVEGEid,'long_name','QVEGE_ratio');
netcdf.putAtt(outid,QVEGEid,'units','unitless');
netcdf.putAtt(outid,QVEGEid,'_FillValue',undef);
netcdf.putAtt(outid,QVEGTid,'long_name','QVEGT_ratio');
netcdf.putAtt(outid,QVEGTid,'units','unitless');
netcdf.putAtt(outid,QVEGTid,'_FillValue',undef);
netcdf.putAtt(outid,QSOILid,'long_name','QSOIL_ratio');
netcdf.putAtt(outid,QSOILid,'units','unitless');
netcdf.putAtt(outid,QSOILid,'_FillValue',undef);
netcdf.close(outid);
Loading...
马建仓 AI 助手
尝试更多
代码解读
代码找茬
代码优化
Matlab
1
https://gitee.com/damowangyang/matlab-scripts-for-gleam.git
git@gitee.com:damowangyang/matlab-scripts-for-gleam.git
damowangyang
matlab-scripts-for-gleam
Matlab scripts for GLEAM
master

搜索帮助