请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
楼主: sysu-lipan6

[程序设计] matlab读取nc文件,并且作等值线降雨图等等

  [复制链接]

新浪微博达人勋

发表于 2015-1-9 18:19:37 | 显示全部楼层

感觉很高大上。。。膜拜
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-1-13 12:44:02 | 显示全部楼层
谢谢分享!!!!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2015-1-13 12:44:07 | 显示全部楼层
谢谢分享!!!!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2015-1-14 11:44:48 | 显示全部楼层
好像很有用的样子 楼主点赞
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-1-14 16:31:22 | 显示全部楼层
谢谢分享,正需要
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-1-14 17:10:54 | 显示全部楼层
% clear everything
clear all;
clc;
path(path,'D:\DataProgram\MatlabLib');
path(path,'D:\DataProgram\MatlabLib\savefigure');

% global 1 by 1 degree data
lon = -179.5:1.0:179.5;
lat = -89.5:1.0:89.5;
n_lon = 360;
n_lat = 180;

% find lon and lat index for amazon, africa, southeast asia, southeast
% america, europe
% (65.5W, 7.5S) Amazon
%( 23.5E,2.5S)  Africa,
% (113.5E,1.5N) Southeast Asia
% (83.5W, 34.5N) Southeast America
% (4.5E, 45.5N) in Europe
lon_amazon = -65.5; lat_amazon = -7.5;
lon_amazon_index = lon == lon_amazon;
lat_amazon_index = lat == lat_amazon;

lon_africa = 23.5; lat_africa = -2.5;
lon_africa_index = lon == lon_africa;
lat_africa_index = lat == lat_africa;

lon_asia = 113.5; lat_asia = 1.5;
lon_asia_index = lon == lon_asia;
lat_asia_index = lat == lat_asia;

lon_america = -83.5; lat_america = 34.5;
lon_america_index = lon == lon_america;
lat_america_index = lat == lat_america;

lon_europe = 4.5; lat_europe = 45.5;
lon_europe_index = lon == lon_europe;
lat_europe_index = lat == lat_europe;

% define input and output directories
inputdir = 'D:\Data\CLM_output\Monthly';
outputdir = 'D:\DataProgram\CLM_output';

% check if input directories exists
if (exist(inputdir,'dir') == 0)
    error('ErrorMsg:Inputdir', ...
          'Input directory %s does not exist.', inputdir);
end;

% check if input directories exists. if not, make one
if (exist(outputdir,'dir') == 0)
    str = sprintf('Output directory %s does not exist.', ...
                  outputdir);
    disp(str);
    disp('The program is making for you.');
    [istatus,mess,messid] = mkdir(outputdir);
    if (istatus == 1)
        disp('The program has made the directory sucessfully');
    else
        error('ErrorMsg:MakeOutputdir', ...
              'The program can not make the directory %s.', outputdir);
    end
end

% ======================================================
% ctl: the CLM3 default scheme
% ======================================================

inputfiles = dir([inputdir, '\', 'CLM3_CTL_60min.clm2.h0.2002*.nc']);
if (length(inputfiles) < 1)
    error('ErrorMsg:InputdirFile', ...
          'There is no required file in the input directory %s.', inputdir);
end
depth = [18 28 45 75 124 204 336 554 913 1137];

% seasonal
for i = 1:length(inputfiles)
    indirfile = [inputdir, '\', inputfiles(i).name];
    if (exist(indirfile,'file') == 2)
        ncid = netcdf.open(indirfile,'NC_NOWRITE');
        if ncid > 0
            
            % rainfall (mm/s) --> (mm/day)
            rainid = netcdf.inqVarID(ncid,'RAIN');
            rain = netcdf.getVar(ncid,rainid);
            rain = double(24.0*3600.0*rain);
            rain(rain > 2.0e+30) = NaN;
            rain_amazon(i) = rain(lon_amazon_index,lat_amazon_index);
            rain_africa(i) = rain(lon_africa_index,lat_africa_index);
            rain_asia(i) = rain(lon_asia_index,lat_asia_index);
            rain_america(i) = rain(lon_america_index,lat_america_index);
            rain_europe(i) = rain(lon_europe_index,lat_europe_index);
            
            
            % incident solar radiation W/m2
            fsdsid = netcdf.inqVarID(ncid,'FSDS');
            fsds = netcdf.getVar(ncid,fsdsid);
            fsds = double(24.0*3600.0*fsds);
            fsds(fsds > 2.0e+30) = NaN;
            fsds_amazon(i) = fsds(lon_amazon_index,lat_amazon_index);
            fsds_africa(i) = fsds(lon_africa_index,lat_africa_index);
            fsds_asia(i) = fsds(lon_asia_index,lat_asia_index);
            fsds_america(i) = fsds(lon_america_index,lat_america_index);
            fsds_europe(i) = fsds(lon_europe_index,lat_europe_index);
            
            % LAI
            elaiid = netcdf.inqVarID(ncid,'ELAI');
            elai = netcdf.getVar(ncid,elaiid);
            elai = double(elai);
            elai(elai > 2.0e+30) = NaN;
            elai_amazon(i) = elai(lon_amazon_index,lat_amazon_index);
            elai_africa(i) = elai(lon_africa_index,lat_africa_index);
            elai_asia(i) = elai(lon_asia_index,lat_asia_index);
            elai_america(i) = elai(lon_america_index,lat_america_index);
            elai_europe(i) = elai(lon_europe_index,lat_europe_index);
            
            % Sensible heat flux (W/m2)
            fshid = netcdf.inqVarID(ncid,'FSH');
            fsh = netcdf.getVar(ncid,fshid);
            fsh = double(fsh);
            fsh(fsh > 2.0e+30) = NaN;
            fsh_amazon_ctl(i) = fsh(lon_amazon_index,lat_amazon_index);
            fsh_africa_ctl(i) = fsh(lon_africa_index,lat_africa_index);
            fsh_asia_ctl(i) = fsh(lon_asia_index,lat_asia_index);
            fsh_america_ctl(i) = fsh(lon_america_index,lat_america_index);
            fsh_europe_ctl(i) = fsh(lon_europe_index,lat_europe_index);
            
            % heat flux for plant transpiration (W/m2)
            fctrid = netcdf.inqVarID(ncid,'FCTR');
            fctr = netcdf.getVar(ncid,fctrid);
            fctr = double(fctr);
            fctr(fctr > 2.0e+30) = NaN;
            
            % heat flux for interception loss (W/m2)
            fcevid = netcdf.inqVarID(ncid,'FCEV');
            fcev = netcdf.getVar(ncid,fcevid);
            fcev = double(fcev);
            fcev(fcev > 2.0e+30) = NaN;            
            
            % heat flux for soil evaporation (W/m2)
            fgevid = netcdf.inqVarID(ncid,'FGEV');
            fgev = netcdf.getVar(ncid,fgevid);
            fgev = double(fgev);
            fgev(fgev > 2.0e+30) = NaN;         
            
            % latent heat flux (W/m2)
            fle = fctr + fcev + fgev;
            fle_amazon_ctl(i) = fle(lon_amazon_index,lat_amazon_index);
            fle_africa_ctl(i) = fle(lon_africa_index,lat_africa_index);
            fle_asia_ctl(i) = fle(lon_asia_index,lat_asia_index);
            fle_america_ctl(i) = fle(lon_america_index,lat_america_index);
            fle_europe_ctl(i) = fle(lon_europe_index,lat_europe_index);
            
            % canopy interception (mm/s) --> (mm/day)
            qvegeid = netcdf.inqVarID(ncid,'QVEGE');
            qvege = netcdf.getVar(ncid,qvegeid);
            qvege = double(24.0*3600.0*qvege);
            qvege(qvege > 2.0e+30) = NaN;
            qvege_amazon_ctl(i) = qvege(lon_amazon_index,lat_amazon_index);
            qvege_africa_ctl(i) = qvege(lon_africa_index,lat_africa_index);
            qvege_asia_ctl(i) = qvege(lon_asia_index,lat_asia_index);
            qvege_america_ctl(i) = qvege(lon_america_index,lat_america_index);
            qvege_europe_ctl(i) = qvege(lon_europe_index,lat_europe_index);           
            
            % plant transpiration (mm/s) --> (mm/day)
            qvegtid = netcdf.inqVarID(ncid,'QVEGT');
            qvegt = netcdf.getVar(ncid,qvegtid);
            qvegt = double(24.0*3600.0*qvegt);
            qvegt(qvegt > 2.0e+30) = NaN;
            qvegt_amazon_ctl(i) = qvegt(lon_amazon_index,lat_amazon_index);
            qvegt_africa_ctl(i) = qvegt(lon_africa_index,lat_africa_index);
            qvegt_asia_ctl(i) = qvegt(lon_asia_index,lat_asia_index);
            qvegt_america_ctl(i) = qvegt(lon_america_index,lat_america_index);
            qvegt_europe_ctl(i) = qvegt(lon_europe_index,lat_europe_index);            
            
            % soil evaporation (mm/s) --> (mm/day)
            qsoilid = netcdf.inqVarID(ncid,'QSOIL');
            qsoil = netcdf.getVar(ncid,qsoilid);
            qsoil = double(24.0*3600.0*qsoil);
            qsoil(qsoil > 2.0e+30) = NaN;
            qsoil_amazon_ctl(i) = qsoil(lon_amazon_index,lat_amazon_index);
            qsoil_africa_ctl(i) = qsoil(lon_africa_index,lat_africa_index);
            qsoil_asia_ctl(i) = qsoil(lon_asia_index,lat_asia_index);
            qsoil_america_ctl(i) = qsoil(lon_america_index,lat_america_index);
            qsoil_europe_ctl(i) = qsoil(lon_europe_index,lat_europe_index);
            
            % total evapotranspiration  (mm/s) --> (mm/day)
            et = qvege + qvegt + qsoil;
            et_amazon_ctl(i) = et(lon_amazon_index,lat_amazon_index);
            et_africa_ctl(i) = et(lon_africa_index,lat_africa_index);
            et_asia_ctl(i) = et(lon_asia_index,lat_asia_index);
            et_america_ctl(i) = et(lon_america_index,lat_america_index);
            et_europe_ctl(i) = et(lon_europe_index,lat_europe_index);            
            
            % surface runoff (mm/s) --> (mm/day)
            qoverid = netcdf.inqVarID(ncid,'QOVER');
            qover = netcdf.getVar(ncid,qoverid);
            qover = double(24.0*3600.0*qover);
            qover(qover > 2.0e+30) = NaN;
            qover_amazon_ctl(i) = qover(lon_amazon_index,lat_amazon_index);
            qover_africa_ctl(i) = qover(lon_africa_index,lat_africa_index);
            qover_asia_ctl(i) = qover(lon_asia_index,lat_asia_index);
            qover_america_ctl(i) = qover(lon_america_index,lat_america_index);
            qover_europe_ctl(i) = qover(lon_europe_index,lat_europe_index);            
            
            % sub-surface runoff (mm/s) --> (mm/day)
            qdraiid = netcdf.inqVarID(ncid,'QDRAI');
            qdrai = netcdf.getVar(ncid,qdraiid);
            qdrai = double(24.0*3600.0*qdrai);
            qdrai(qdrai > 2.0e+30) = NaN;
            qdrai_amazon_ctl(i) = qdrai(lon_amazon_index,lat_amazon_index);
            qdrai_africa_ctl(i) = qdrai(lon_africa_index,lat_africa_index);
            qdrai_asia_ctl(i) = qdrai(lon_asia_index,lat_asia_index);
            qdrai_america_ctl(i) = qdrai(lon_america_index,lat_america_index);
            qdrai_europe_ctl(i) = qdrai(lon_europe_index,lat_europe_index);            
            
            % runoff (mm/s) --> (mm/day)
            runoff = qover + qdrai;
            runoff_amazon_ctl(i) = runoff(lon_amazon_index,lat_amazon_index);
            runoff_africa_ctl(i) = runoff(lon_africa_index,lat_africa_index);
            runoff_asia_ctl(i) = runoff(lon_asia_index,lat_asia_index);
            runoff_america_ctl(i) = runoff(lon_america_index,lat_america_index);
            runoff_europe_ctl(i) = runoff(lon_europe_index,lat_europe_index);            
            
            % ground temperature (K)
            tgid = netcdf.inqVarID(ncid,'TG');
            tg = netcdf.getVar(ncid,tgid);
            tg = double(tg);
            tg(tg > 2.0e+30) = NaN;
            tg_amazon_ctl(i) = tg(lon_amazon_index,lat_amazon_index);
            tg_africa_ctl(i) = tg(lon_africa_index,lat_africa_index);
            tg_asia_ctl(i) = tg(lon_asia_index,lat_asia_index);
            tg_america_ctl(i) = tg(lon_america_index,lat_america_index);
            tg_europe_ctl(i) = tg(lon_europe_index,lat_europe_index);            
            
            % vegetation temperature (K)
            tvid = netcdf.inqVarID(ncid,'TV');
            tv = netcdf.getVar(ncid,tvid);
            tv = double(tv);
            tv(tv > 2.0e+30) = NaN;
            tv_amazon_ctl(i) = tv(lon_amazon_index,lat_amazon_index);
            tv_africa_ctl(i) = tv(lon_africa_index,lat_africa_index);
            tv_asia_ctl(i) = tv(lon_asia_index,lat_asia_index);
            tv_america_ctl(i) = tv(lon_america_index,lat_america_index);
            tv_europe_ctl(i) = tv(lon_europe_index,lat_europe_index);
            
            % soil water (kg/m2) == mm
            soilliqid = netcdf.inqVarID(ncid,'SOILLIQ');
            soilliq = netcdf.getVar(ncid,soilliqid);
            soilliq = double(soilliq);
            soilliq(soilliq > 2.0e+30) = NaN;
            
            soiliceid = netcdf.inqVarID(ncid,'SOILICE');
            soilice = netcdf.getVar(ncid,soiliceid);
            soilice = double(soilice);
            soilice(soilice > 2.0e+30) = NaN;
            
            soil = soilliq + soilice;
            soil(soil > 2.0e+30) = NaN;
            
            % total soil moisture for 10 layers (3.4m)
            soil_total = sum(soil,3,'native');
            soil_amazon_ctl(i) = soil_total(lon_amazon_index,lat_amazon_index);
            soil_africa_ctl(i) = soil_total(lon_africa_index,lat_africa_index);
            soil_asia_ctl(i) = soil_total(lon_asia_index,lat_asia_index);
            soil_america_ctl(i) = soil_total(lon_america_index,lat_america_index);
            soil_europe_ctl(i) = soil_total(lon_europe_index,lat_europe_index);
            
            soil_1m_total = sum(soil(:,:,1:7),3,'native') + soil(:,:,8)*0.1711;
            soil_1m_amazon_ctl(i) = soil_1m_total(lon_amazon_index,lat_amazon_index);
            soil_1m_africa_ctl(i) = soil_1m_total(lon_africa_index,lat_africa_index);
            soil_1m_asia_ctl(i) = soil_1m_total(lon_asia_index,lat_asia_index);
            soil_1m_america_ctl(i) = soil_1m_total(lon_america_index,lat_america_index);
            soil_1m_europe_ctl(i) = soil_1m_total(lon_europe_index,lat_europe_index);
            
            % soil temperature (K)
            tsoiid = netcdf.inqVarID(ncid,'TSOI');
            tsoi = netcdf.getVar(ncid,tsoiid);
            tsoi = double(tsoi);
            tsoi(tsoi > 2.0e+30) = NaN;
            
            tsoi_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
                      + tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
                      + tsoi(:,:,7)*336 + tsoi(:,:,8)*554 + tsoi(:,:,9)*913 + tsoi(:,:,9)*1137)/3400;
            tsoi_amazon_ctl(i) = tsoi_ave(lon_amazon_index,lat_amazon_index);
            tsoi_africa_ctl(i) = tsoi_ave(lon_africa_index,lat_africa_index);
            tsoi_asia_ctl(i) = tsoi_ave(lon_asia_index,lat_asia_index);
            tsoi_america_ctl(i) = tsoi_ave(lon_america_index,lat_america_index);
            tsoi_europe_ctl(i) = tsoi_ave(lon_europe_index,lat_europe_index);                 
                  
            tsoi_1m_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
                      + tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
                      + tsoi(:,:,7)*336 + tsoi(:,:,8)*171)/1000;
            tsoi_1m_amazon_ctl(i) = tsoi_1m_ave(lon_amazon_index,lat_amazon_index);
            tsoi_1m_africa_ctl(i) = tsoi_1m_ave(lon_africa_index,lat_africa_index);
            tsoi_1m_asia_ctl(i) = tsoi_1m_ave(lon_asia_index,lat_asia_index);
            tsoi_1m_america_ctl(i) = tsoi_1m_ave(lon_america_index,lat_america_index);
            tsoi_1m_europe_ctl(i) = tsoi_1m_ave(lon_europe_index,lat_europe_index);            
            
            netcdf.close(ncid);
        else
            error('ErrorMsg:indirfileOpen', ...
                  'File %s can not be opened.', indirfile);
        end
        
    else
        error('ErrorMsg:indirfileExistOrNot', ...
              'File %s does not exist.', indirfile);
    end
end

% ======================================================
% v2: the improved scheme
% ======================================================

inputfiles = dir([inputdir, '\', 'CLM3_EXP3_v2_60min.clm2.h0.2002*.nc']);
if (length(inputfiles) < 1)
    error('ErrorMsg:InputdirFile', ...
          'There is no required file in the input directory %s.', inputdir);
end

% seasonal
for i = 1:length(inputfiles)
    indirfile = [inputdir, '\', inputfiles(i).name];
    if (exist(indirfile,'file') == 2)
        ncid = netcdf.open(indirfile,'NC_NOWRITE');
        if ncid > 0
            
            % rainfall (mm/s) --> (mm/day)
%             rainid = netcdf.inqVarID(ncid,'RAIN');
%             rain = netcdf.getVar(ncid,rainid);
%             rain = double(24.0*3600.0*rain);
%             rain(rain > 2.0e+30) = NaN;
%             rain_amazon(i) = rain(lon_amazon_index,lat_amazon_index);
%             rain_africa(i) = rain(lon_africa_index,lat_africa_index);
%             rain_asia(i) = rain(lon_asia_index,lat_asia_index);
%             rain_america(i) = rain(lon_america_index,lat_america_index);
%             rain_europe(i) = rain(lon_europe_index,lat_europe_index);
%            
%             % LAI
%             elaiid = netcdf.inqVarID(ncid,'ELAI');
%             elai = netcdf.getVar(ncid,elaiid);
%             elai = double(elai);
%             elai(elai > 2.0e+30) = NaN;
%             elai_amazon(i) = elai(lon_amazon_index,lat_amazon_index);
%             elai_africa(i) = elai(lon_africa_index,lat_africa_index);
%             elai_asia(i) = elai(lon_asia_index,lat_asia_index);
%             elai_america(i) = elai(lon_america_index,lat_america_index);
%             elai_europe(i) = elai(lon_europe_index,lat_europe_index);
            
            % Sensible heat flux (W/m2)
            fshid = netcdf.inqVarID(ncid,'FSH');
            fsh = netcdf.getVar(ncid,fshid);
            fsh = double(fsh);
            fsh(fsh > 2.0e+30) = NaN;
            fsh_amazon_exp(i) = fsh(lon_amazon_index,lat_amazon_index);
            fsh_africa_exp(i) = fsh(lon_africa_index,lat_africa_index);
            fsh_asia_exp(i) = fsh(lon_asia_index,lat_asia_index);
            fsh_america_exp(i) = fsh(lon_america_index,lat_america_index);
            fsh_europe_exp(i) = fsh(lon_europe_index,lat_europe_index);
            
            % heat flux for plant transpiration (W/m2)
            fctrid = netcdf.inqVarID(ncid,'FCTR');
            fctr = netcdf.getVar(ncid,fctrid);
            fctr = double(fctr);
            fctr(fctr > 2.0e+30) = NaN;
            
            % heat flux for interception loss (W/m2)
            fcevid = netcdf.inqVarID(ncid,'FCEV');
            fcev = netcdf.getVar(ncid,fcevid);
            fcev = double(fcev);
            fcev(fcev > 2.0e+30) = NaN;            
            
            % heat flux for soil evaporation (W/m2)
            fgevid = netcdf.inqVarID(ncid,'FGEV');
            fgev = netcdf.getVar(ncid,fgevid);
            fgev = double(fgev);
            fgev(fgev > 2.0e+30) = NaN;         
            
            % latent heat flux (W/m2)
            fle = fctr + fcev + fgev;
            fle_amazon_exp(i) = fle(lon_amazon_index,lat_amazon_index);
            fle_africa_exp(i) = fle(lon_africa_index,lat_africa_index);
            fle_asia_exp(i) = fle(lon_asia_index,lat_asia_index);
            fle_america_exp(i) = fle(lon_america_index,lat_america_index);
            fle_europe_exp(i) = fle(lon_europe_index,lat_europe_index);
            
            % canopy interception (mm/s) --> (mm/day)
            qvegeid = netcdf.inqVarID(ncid,'QVEGE');
            qvege = netcdf.getVar(ncid,qvegeid);
            qvege = double(24.0*3600.0*qvege);
            qvege(qvege > 2.0e+30) = NaN;
            qvege_amazon_exp(i) = qvege(lon_amazon_index,lat_amazon_index);
            qvege_africa_exp(i) = qvege(lon_africa_index,lat_africa_index);
            qvege_asia_exp(i) = qvege(lon_asia_index,lat_asia_index);
            qvege_america_exp(i) = qvege(lon_america_index,lat_america_index);
            qvege_europe_exp(i) = qvege(lon_europe_index,lat_europe_index);           
            
            % plant transpiration (mm/s) --> (mm/day)
            qvegtid = netcdf.inqVarID(ncid,'QVEGT');
            qvegt = netcdf.getVar(ncid,qvegtid);
            qvegt = double(24.0*3600.0*qvegt);
            qvegt(qvegt > 2.0e+30) = NaN;
            qvegt_amazon_exp(i) = qvegt(lon_amazon_index,lat_amazon_index);
            qvegt_africa_exp(i) = qvegt(lon_africa_index,lat_africa_index);
            qvegt_asia_exp(i) = qvegt(lon_asia_index,lat_asia_index);
            qvegt_america_exp(i) = qvegt(lon_america_index,lat_america_index);
            qvegt_europe_exp(i) = qvegt(lon_europe_index,lat_europe_index);            
            
            % soil evaporation (mm/s) --> (mm/day)
            qsoilid = netcdf.inqVarID(ncid,'QSOIL');
            qsoil = netcdf.getVar(ncid,qsoilid);
            qsoil = double(24.0*3600.0*qsoil);
            qsoil(qsoil > 2.0e+30) = NaN;
            qsoil_amazon_exp(i) = qsoil(lon_amazon_index,lat_amazon_index);
            qsoil_africa_exp(i) = qsoil(lon_africa_index,lat_africa_index);
            qsoil_asia_exp(i) = qsoil(lon_asia_index,lat_asia_index);
            qsoil_america_exp(i) = qsoil(lon_america_index,lat_america_index);
            qsoil_europe_exp(i) = qsoil(lon_europe_index,lat_europe_index);
            
            % total evapotranspiration  (mm/s) --> (mm/day)
            et = qvege + qvegt + qsoil;
            et_amazon_exp(i) = et(lon_amazon_index,lat_amazon_index);
            et_africa_exp(i) = et(lon_africa_index,lat_africa_index);
            et_asia_exp(i) = et(lon_asia_index,lat_asia_index);
            et_america_exp(i) = et(lon_america_index,lat_america_index);
            et_europe_exp(i) = et(lon_europe_index,lat_europe_index);            
            
            % surface runoff (mm/s) --> (mm/day)
            qoverid = netcdf.inqVarID(ncid,'QOVER');
            qover = netcdf.getVar(ncid,qoverid);
            qover = double(24.0*3600.0*qover);
            qover(qover > 2.0e+30) = NaN;
            qover_amazon_exp(i) = qover(lon_amazon_index,lat_amazon_index);
            qover_africa_exp(i) = qover(lon_africa_index,lat_africa_index);
            qover_asia_exp(i) = qover(lon_asia_index,lat_asia_index);
            qover_america_exp(i) = qover(lon_america_index,lat_america_index);
            qover_europe_exp(i) = qover(lon_europe_index,lat_europe_index);            
            
            % sub-surface runoff (mm/s) --> (mm/day)
            qdraiid = netcdf.inqVarID(ncid,'QDRAI');
            qdrai = netcdf.getVar(ncid,qdraiid);
            qdrai = double(24.0*3600.0*qdrai);
            qdrai(qdrai > 2.0e+30) = NaN;
            qdrai_amazon_exp(i) = qdrai(lon_amazon_index,lat_amazon_index);
            qdrai_africa_exp(i) = qdrai(lon_africa_index,lat_africa_index);
            qdrai_asia_exp(i) = qdrai(lon_asia_index,lat_asia_index);
            qdrai_america_exp(i) = qdrai(lon_america_index,lat_america_index);
            qdrai_europe_exp(i) = qdrai(lon_europe_index,lat_europe_index);            
            
            % runoff (mm/s) --> (mm/day)
            runoff = qover + qdrai;
            runoff_amazon_exp(i) = runoff(lon_amazon_index,lat_amazon_index);
            runoff_africa_exp(i) = runoff(lon_africa_index,lat_africa_index);
            runoff_asia_exp(i) = runoff(lon_asia_index,lat_asia_index);
            runoff_america_exp(i) = runoff(lon_america_index,lat_america_index);
            runoff_europe_exp(i) = runoff(lon_europe_index,lat_europe_index);            
            
            % ground temperature (K)
            tgid = netcdf.inqVarID(ncid,'TG');
            tg = netcdf.getVar(ncid,tgid);
            tg = double(tg);
            tg(tg > 2.0e+30) = NaN;
            tg_amazon_exp(i) = tg(lon_amazon_index,lat_amazon_index);
            tg_africa_exp(i) = tg(lon_africa_index,lat_africa_index);
            tg_asia_exp(i) = tg(lon_asia_index,lat_asia_index);
            tg_america_exp(i) = tg(lon_america_index,lat_america_index);
            tg_europe_exp(i) = tg(lon_europe_index,lat_europe_index);            
            
            % vegetation temperature (K)
            tvid = netcdf.inqVarID(ncid,'TV');
            tv = netcdf.getVar(ncid,tvid);
            tv = double(tv);
            tv(tv > 2.0e+30) = NaN;
            tv_amazon_exp(i) = tv(lon_amazon_index,lat_amazon_index);
            tv_africa_exp(i) = tv(lon_africa_index,lat_africa_index);
            tv_asia_exp(i) = tv(lon_asia_index,lat_asia_index);
            tv_america_exp(i) = tv(lon_america_index,lat_america_index);
            tv_europe_exp(i) = tv(lon_europe_index,lat_europe_index);
            
            % soil water (kg/m2) == mm
            soilliqid = netcdf.inqVarID(ncid,'SOILLIQ');
            soilliq = netcdf.getVar(ncid,soilliqid);
            soilliq = double(soilliq);
            soilliq(soilliq > 2.0e+30) = NaN;
            
            soiliceid = netcdf.inqVarID(ncid,'SOILICE');
            soilice = netcdf.getVar(ncid,soiliceid);
            soilice = double(soilice);
            soilice(soilice > 2.0e+30) = NaN;
            
            soil = soilliq + soilice;
            soil(soil > 2.0e+30) = NaN;
            
            % total soil moisture for 10 layers (3.4m)
            soil_total = sum(soil,3,'native');
            soil_amazon_exp(i) = soil_total(lon_amazon_index,lat_amazon_index);
            soil_africa_exp(i) = soil_total(lon_africa_index,lat_africa_index);
            soil_asia_exp(i) = soil_total(lon_asia_index,lat_asia_index);
            soil_america_exp(i) = soil_total(lon_america_index,lat_america_index);
            soil_europe_exp(i) = soil_total(lon_europe_index,lat_europe_index);
            
            soil_1m_total = sum(soil(:,:,1:7),3,'native') + soil(:,:,8)*0.1711;
            soil_1m_amazon_exp(i) = soil_1m_total(lon_amazon_index,lat_amazon_index);
            soil_1m_africa_exp(i) = soil_1m_total(lon_africa_index,lat_africa_index);
            soil_1m_asia_exp(i) = soil_1m_total(lon_asia_index,lat_asia_index);
            soil_1m_america_exp(i) = soil_1m_total(lon_america_index,lat_america_index);
            soil_1m_europe_exp(i) = soil_1m_total(lon_europe_index,lat_europe_index);
            
            % soil temperature (K)
            tsoiid = netcdf.inqVarID(ncid,'TSOI');
            tsoi = netcdf.getVar(ncid,tsoiid);
            tsoi = double(tsoi);
            tsoi(tsoi > 2.0e+30) = NaN;
            
            tsoi_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
                      + tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
                      + tsoi(:,:,7)*336 + tsoi(:,:,8)*554 + tsoi(:,:,9)*913 + tsoi(:,:,9)*1137)/3400;
            tsoi_amazon_exp(i) = tsoi_ave(lon_amazon_index,lat_amazon_index);
            tsoi_africa_exp(i) = tsoi_ave(lon_africa_index,lat_africa_index);
            tsoi_asia_exp(i) = tsoi_ave(lon_asia_index,lat_asia_index);
            tsoi_america_exp(i) = tsoi_ave(lon_america_index,lat_america_index);
            tsoi_europe_exp(i) = tsoi_ave(lon_europe_index,lat_europe_index);                 
                  
            tsoi_1m_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
                      + tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
                      + tsoi(:,:,7)*336 + tsoi(:,:,8)*171)/1000;
            tsoi_1m_amazon_exp(i) = tsoi_1m_ave(lon_amazon_index,lat_amazon_index);
            tsoi_1m_africa_exp(i) = tsoi_1m_ave(lon_africa_index,lat_africa_index);
            tsoi_1m_asia_exp(i) = tsoi_1m_ave(lon_asia_index,lat_asia_index);
            tsoi_1m_america_exp(i) = tsoi_1m_ave(lon_america_index,lat_america_index);
            tsoi_1m_europe_exp(i) = tsoi_1m_ave(lon_europe_index,lat_europe_index);            
            
            netcdf.close(ncid);
        else
            error('ErrorMsg:indirfileOpen', ...
                  'File %s can not be opened.', indirfile);
        end
        
    else
        error('ErrorMsg:indirfileExistOrNot', ...
              'File %s does not exist.', indirfile);
    end
end

% rain coverage
inputdir = 'D:\Data\CLM_output\ForcingAndObservedRain';
inputfiles = dir([inputdir, '\', 'raincoverage-*.txt']);
if (length(inputfiles) < 1)
    error('ErrorMsg:InputdirFile', ...
          'There is no required file in the input directory %s.', inputdir);
end

for i = 1:length(inputfiles)
    indirfile = [inputdir, '\', inputfiles(i).name];
    fid = fopen(indirfile,'r');
    raincoverage = fscanf(fid, '%f', [360, 180]);
    raincoverage_amazon(i) = raincoverage(lon_amazon_index,lat_amazon_index);
    raincoverage_africa(i) = raincoverage(lon_africa_index,lat_africa_index);
    raincoverage_asia(i) = raincoverage(lon_asia_index,lat_asia_index);
    raincoverage_america(i) = raincoverage(lon_america_index,lat_america_index);
    raincoverage_europe(i) = raincoverage(lon_europe_index,lat_europe_index);
    fclose(fid);
end


% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% &&&&&&&&&&&&&&&& plot &&&&&&&&&&&&&&&&&&&&&&&&&&&&
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

xx = 0:1:11;
scrsz = get(0,'ScreenSize');
figure('Position',[5 55 scrsz(3)-10 scrsz(4)-115]);
subplot1(2,5,'Gap',[0.04 0.07],'XTickL','Margin','YTickL','All','FontS',8,'Min',[0.03 0.1],'Max',[0.99 0.99]);

subplot1(1);
plot(xx,rain_amazon,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
ylabel('Rainfall (mm/day)','FontSize',9);
title('AM-65.5W-7.5S','fontsize',14);
text(0.5,11.4,'a1','FontSize',7,'FontWeight','bold');
grid on;

subplot1(2);
plot(xx,rain_africa,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('AF-23.5E-2.5S','fontsize',14);
text(0.5,7.5,'a2','FontSize',7,'FontWeight','bold');
grid on;

subplot1(3);
plot(xx,rain_asia,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('MA-113.5E-1.5N','fontsize',14);
text(0.5,9.5,'a3','FontSize',7,'FontWeight','bold');
grid on;

subplot1(4);
plot(xx,rain_america,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('US-83.5W-34.5N','fontsize',14);
text(0.5,14,'a4','FontSize',7,'FontWeight','bold');
grid on;

subplot1(5);
plot(xx,rain_europe,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('EU-4.5E-45.5N','fontsize',14);
text(0.5,5.5,'a5','FontSize',7,'FontWeight','bold');
grid on;

subplot1(6);
plot(xx,elai_amazon,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:12);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
ylabel('Leaf area index','FontSize',9);
text(0.5,5.9,'b1','FontSize',7,'FontWeight','bold');
grid on;

subplot1(7);
plot(xx,elai_africa,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:12);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,5.9,'b2','FontSize',7,'FontWeight','bold');
grid on;

subplot1(8);
plot(xx,elai_asia,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,5.35,'b3','FontSize',7,'FontWeight','bold');
grid on;

subplot1(9);
plot(xx,elai_america,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,4.75,'b4','FontSize',7,'FontWeight','bold');
grid on;

subplot1(10);
plot(xx,elai_europe,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,2.93,'b5','FontSize',7,'FontWeight','bold');
grid on;

set(gcf,'PaperPositionMode','auto')
savefigure('Revised_Figure3_Seasonal_rain_lai.tif','PrintOptions','-dtiff', 'Size',[6.8 3],'LineWidth',1)
disp('done');
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-1-14 17:12:03 | 显示全部楼层
第三份代码图我找不到了……惭愧,但是也有很高是艺术价值,具体说来,有几层楼那么高
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-1-15 13:08:24 | 显示全部楼层
楼主方便告诉您的数据是哪里下的吗?谢谢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-1-16 14:15:13 | 显示全部楼层
woniu123 发表于 2015-1-15 13:08
楼主方便告诉您的数据是哪里下的吗?谢谢

老师给的,他去nc那个官网下的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-1-16 21:11:20 | 显示全部楼层
能介绍下是啥程序不?
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表