matlab读取nc文件，并且作等值线降雨图等等

sysu-lipan6

精品啊，老师给的参考程序，冒死上传，被发现要出事，后续应该还有

wuzzh_gd

那我瞅瞅。。

更大华

有这么严重啊

893116547

这么严重啊，那你还传，

sysu-lipan6

% clear everything
clear all;
clc;
path(path,'D:\DataProgram\MatlabLib');
path(path,'D:\DataProgram\MatlabLib\savefigure');

% define longitude and lantitude
% the spatial resolution is 1 by 1 degree globally
x = -179.5:1.0:179.5;
y = -89.5:1.0:89.5;
n_lat = 180;
n_lon = 360;

% 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.*.nc']);
if (length(inputfiles) < 1)
    error('ErrorMsg:InputdirFile', ...
          'There is no required file in the input directory %s.', inputdir);
end

% to speed up running, define size of array
lenthfiles = zeros(n_lon,n_lat,'double');
totl_et = zeros(n_lon,n_lat,'double');
et_ctl = zeros(n_lon,n_lat,'double');

% calculate total rain, snow, and canopy evaporation
for i = 1:length(inputfiles)
    indirfile = [inputdir, '/', inputfiles(i).name];
    if (exist(indirfile,'file') == 2)
        ncid = netcdf.open(indirfile,'NC_NOWRITE');   
        if ncid > 0     
            snowid = netcdf.inqVarID(ncid,'SNOW');
            snow = netcdf.getVar(ncid,snowid);
            snow = double(24.0*3600.0*snow);
            snow(snow > 2.0e+30) = NaN;
            
            rainid = netcdf.inqVarID(ncid,'RAIN');
            rain = netcdf.getVar(ncid,rainid);
            rain = double(24.0*3600.0*rain);
            rain(rain > 2.0e+30) = NaN;
            
            % canopy interception
            qvegeid = netcdf.inqVarID(ncid,'QVEGE');
            qvege = netcdf.getVar(ncid,qvegeid);
            qvege = double(24.0*3600.0*qvege);
            qvege(qvege > 2.0e+30) = NaN;
            
            % plant transpiration
            qvegtid = netcdf.inqVarID(ncid,'QVEGT');
            qvegt = netcdf.getVar(ncid,qvegtid);
            qvegt = double(24.0*3600.0*qvegt);
            qvegt(qvegt > 2.0e+30) = NaN;
            
            % soil evaporation
            qsoilid = netcdf.inqVarID(ncid,'QSOIL');
            qsoil = netcdf.getVar(ncid,qsoilid);
            qsoil = double(24.0*3600.0*qsoil);
            qsoil(qsoil > 2.0e+30) = NaN;
            
            et = qvege + qvegt + qsoil;
            totl_et = totl_et + et;
            lenthfiles = lenthfiles +1.0;
        else
            error('ErrorMsg:indirfileOpen', ...
                  'File %s can not be opened.', indirfile);
        end
        netcdf.close(ncid);
    else
        error('ErrorMsg:indirfileExistOrNot', ...
              'File %s does not exist.', indirfile);
    end
end
      
% calculate averages of rain, canopy evaporation
% and ratio of canopy interception to rain
for i_lat = 1:n_lat
    for i_lon = 1:n_lon
        if (lenthfiles(i_lon,i_lat) > 0.)
            et_ctl(i_lon,i_lat) = totl_et(i_lon,i_lat)/lenthfiles(i_lon,i_lat);
        else
            et_ctl(i_lon,i_lat) = NaN;
        end
    end
end


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

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

% to speed up running, define size of array
lenthfiles = zeros(n_lon,n_lat,'double');
totl_et = zeros(n_lon,n_lat,'double');
et_exp = zeros(n_lon,n_lat,'double');


% calculate total rain, snow, and canopy evaporation
for i = 1:length(inputfiles)
    indirfile = [inputdir, '/', inputfiles(i).name];
    if (exist(indirfile,'file') == 2)
        ncid = netcdf.open(indirfile,'NC_NOWRITE');   
        if ncid > 0     
            snowid = netcdf.inqVarID(ncid,'SNOW');
            snow = netcdf.getVar(ncid,snowid);
            snow = double(24.0*3600.0*snow);
            snow(snow > 2.0e+30) = NaN;
            
            rainid = netcdf.inqVarID(ncid,'RAIN');
            rain = netcdf.getVar(ncid,rainid);
            rain = double(24.0*3600.0*rain);
            rain(rain > 2.0e+30) = NaN;
            
            % canopy interception
            qvegeid = netcdf.inqVarID(ncid,'QVEGE');
            qvege = netcdf.getVar(ncid,qvegeid);
            qvege = double(24.0*3600.0*qvege);
            qvege(qvege > 2.0e+30) = NaN;
            
            % plant transpiration
            qvegtid = netcdf.inqVarID(ncid,'QVEGT');
            qvegt = netcdf.getVar(ncid,qvegtid);
            qvegt = double(24.0*3600.0*qvegt);
            qvegt(qvegt > 2.0e+30) = NaN;
            
            % soil evaporation
            qsoilid = netcdf.inqVarID(ncid,'QSOIL');
            qsoil = netcdf.getVar(ncid,qsoilid);
            qsoil = double(24.0*3600.0*qsoil);
            qsoil(qsoil > 2.0e+30) = NaN;
            
            et = qvege + qvegt + qsoil;
            totl_et = totl_et + et;
            lenthfiles = lenthfiles +1.0;
        else
            error('ErrorMsg:indirfileOpen', ...
                  'File %s can not be opened.', indirfile);
        end
        netcdf.close(ncid);
    else
        error('ErrorMsg:indirfileExistOrNot', ...
              'File %s does not exist.', indirfile);
    end
end
   
% calculate averages of rain, canopy evaporation
% and ratio of canopy interception to rain
for i_lat = 1:n_lat
    for i_lon = 1:n_lon
        if (lenthfiles(i_lon,i_lat) > 0.)
            et_exp(i_lon,i_lat) = totl_et(i_lon,i_lat)/lenthfiles(i_lon,i_lat);
        else
            et_exp(i_lon,i_lat) = NaN;
        end
    end
end

% difference between exp and ctl
et_exp_ctl = et_exp - et_ctl;

% keep range -90-60, 80-90
et_ctl(:,1:30) = nan; et_ctl(:,171:180) = nan;
et_exp(:,1:30) = nan; et_exp(:,171:180) = nan;
et_exp_ctl(:,1:30) = nan; et_exp_ctl(:,171:180) = nan;

% ======================================================
% MODIS
% ======================================================
inputdir = 'D:\Data\GlobalET\MOD16A2_MONTHLY.MERRA_GMAO_1kmALB';
outputdir = 'D:\DataProgram\CLM_output';
nRowMODIS = 280; nColMODIS = 720;
lat1MODIS = -60; lat2MODIS = 80; lon1MODIS = -180; lon2MODIS = 180;
yStart = 2000; yEnd = 2005;

% integrate to 1 year
AVEETMODIS = zeros(nRowMODIS,nColMODIS);
for iYear = yStart:yEnd
    for iMon = 1:12
        iMonStr = num2str(iMon);
        if iMon < 10, iMonStr = ['0',num2str(iMon)];end
        yearMonStr = [num2str(iYear), 'M', iMonStr];
        inFile = ['MOD16A2_ET_0.5deg_GEO_', yearMonStr, '.tif'];
        inDirFile = [inputdir, '\', inFile];
        [ETMODIS, cmap, refmat, bbox] = geotiffread(inDirFile);
        ETMODIS(ETMODIS>30000) = nan;
        AVEETMODIS = AVEETMODIS + double(ETMODIS);
    end
end
nYear = (yEnd - yStart + 1);

% convert unit from 0.1mm to 1mm
AVEETMODIS = AVEETMODIS/nYear/365.25/10.0;

% integrate to 1 degree
AVEETMODIS1D = zeros(nRowMODIS/2,nColMODIS/2);
for i = 1:nRowMODIS/2   
    for j = 1:nColMODIS/2
        iStart = (i-1)*2 + 1;
        iEnd = (i-1)*2 + 2;
        jStart = (j-1)*2 + 1;
        jEnd = (j-1)*2 + 2;
        ETMODIS1D = 0;
        nETMODIS1D = 0;
        for ii = iStart:iEnd
            for jj = jStart:jEnd
                ETMODIS1D = ETMODIS1D + AVEETMODIS(ii,jj);
                nETMODIS1D =  nETMODIS1D + 1;
            end
        end
        AVEETMODIS1D(i,j) = ETMODIS1D/nETMODIS1D;
    end
end

% put small -60-80 range into big (-90-90) map
AVEETMODIS1DGlabal = nan(360,180);
AVEETMODIS1DGlabal(:,31:170) = fliplr(AVEETMODIS1D');
AVEETMODIS1DGlabal(isnan(et_ctl)) = nan;

et_ctl_modis = et_ctl - AVEETMODIS1DGlabal;
et_exp_modis = et_exp - AVEETMODIS1DGlabal;

% regional average
ave_et_ctl = mean(mean(et_ctl(isnan(et_ctl)==0)));
ave_et_exp = mean(mean(et_exp(isnan(et_exp)==0)));
ave_AVEETMODIS1DGlabal = mean(mean(AVEETMODIS1DGlabal(isnan(AVEETMODIS1DGlabal)==0)));
ave_et_exp_ctl = ave_et_exp - ave_et_ctl;
ave_et_ctl_modis = ave_et_ctl - ave_AVEETMODIS1DGlabal;
ave_et_exp_modis = ave_et_exp - ave_AVEETMODIS1DGlabal;

% ======================================================
% plot
% ======================================================
scrsz = get(0,'ScreenSize');
figure_h = figure('Position',[5 55 scrsz(3)-10 scrsz(4)-115]);
% rgbc = [127 0 255; 0 0 255; 0 128 255; 0 255 255; 0 255 128; 128 255 0; 255 255 0; 255 128 0; 255 0 0; 153 0 0;...
%        0 0 200; 55 55 255; 110 110 255; 165 165 255; 220 220 255; 255 220 220; 255 165 165; 255 110 110; 255 55 55; 200 0 0]/255;
rgbc = [153 0 0; 255 0 0; 255 128 0; 255 255 0; 128 255 0; 0 255 128; 0 255 255; 0 128 255; 0 0 255; 127 0 255; ...
        200 0 0; 255 55 55; 255 110 110; 255 165 165; 255 220 220; 220 220 255; 165 165 255; 110 110 255; 55 55 255; 0 0 200]/255;
colormap(rgbc);
subplot1(2,2,'Gap',[0.05 0.08],'XTickL','All','YTickL','All','FontS',6,'Min',[0.03 0.05],'Max',[0.99 0.99]);

% colorbar level
cMin = 0;
cMax = 6;
cTick = (cMax - cMin)/10;
cxisMin = cMin;
cxisMax = cMax + (cMax - cMin);
for iLabel = 1:9
    xLabelNum = cMin + iLabel*cTick;
    xLabelStr{iLabel} = num2str(xLabelNum);
end

% ctl
subplot1(1)
v = cMin:cTick:cMax;
et_ctl(et_ctl > cMax) = cMax;
[~, h] = contourf(x,y,et_ctl',v);hold on;
set(h,'edgecolor','none');
axis equal;
caxis([cxisMin cxisMax]);
h = colorbar('horiz');
set(h,'XLim',[cMin cMax]);
set(h,'XTick',cMin+cTick:cTick:cMax-cTick);
set(h,'XTickLabel',{xLabelStr{1},xLabelStr{2},xLabelStr{3},xLabelStr{4},...
    xLabelStr{5},xLabelStr{6},xLabelStr{7},xLabelStr{8},xLabelStr{9}});
set(h,'position', [0.03 0.52 0.43 0.015]);
set(h,'fontsize',11);
text(-173,83,'a','FontSize',7,'FontWeight','bold');
load coast;
plot(long,lat,'k-','LineWidth',0.6); hold off;
set(gca,'XLim',[-180 180],'YLim',[-90 90]);
set(gca,'XTick',-120:60:120);
set(gca,'XTickLabel',{'','','','',''});
set(gca,'YTick',-60:30:60);
set(gca,'YTickLabel',{'60S','30S','EQ','30N','60N'});
titl = ['CTL ',num2str2dig(ave_et_ctl)];
title(titl,'fontsize',14);
grid on;

% colorbar level
cMin = -2.5;
cMax = 2.5;
cTick = (cMax - cMin)/10;
cxisMin = cMin - (cMax - cMin);
cxisMax = cMax;
for iLabel = 1:9
    xLabelNum = cMin + iLabel*cTick;
    xLabelStr{iLabel} = num2str(xLabelNum);
end

% exp-ctl
subplot1(2)
v = cMin:cTick:cMax;
et_exp_ctl(et_exp_ctl > cMax) = cMax;
[~, h] = contourf(x,y,et_exp_ctl',v);hold on;
set(h,'edgecolor','none');
axis equal;
caxis([cxisMin cxisMax]);
h = colorbar('horiz');
set(h,'XLim',[cMin cMax]);
set(h,'XTick',cMin+cTick:cTick:cMax-cTick);
set(h,'XTickLabel',{xLabelStr{1},xLabelStr{2},xLabelStr{3},xLabelStr{4},...
    xLabelStr{5},xLabelStr{6},xLabelStr{7},xLabelStr{8},xLabelStr{9}});
set(h,'position', [0.51 0.52 0.43 0.015]);
set(h,'fontsize',11);
text(-173,83,'b','FontSize',7,'FontWeight','bold');
load coast;
plot(long,lat,'k-','LineWidth',0.6); hold off;
set(gca,'XLim',[-180 180],'YLim',[-90 90]);
set(gca,'XTick',-120:60:120);
set(gca,'XTickLabel',{'','','','',''});
set(gca,'YTick',-60:30:60);
set(gca,'YTickLabel',{'','','','',''});
titl = ['EXP—CTL ',num2str2dig(ave_et_exp_ctl)];
title(titl,'fontsize',14);
grid on;

% ctl-modis
subplot1(3)
v = cMin:cTick:cMax;
et_ctl_modis(et_ctl_modis > cMax) = cMax;
[~, h] = contourf(x,y,et_ctl_modis',v);hold on;
set(h,'edgecolor','none');
axis equal;
caxis([cxisMin cxisMax]);
h = colorbar('horiz');
set(h,'XLim',[cMin cMax]);
set(h,'XTick',cMin+cTick:cTick:cMax-cTick);
set(h,'XTickLabel',{xLabelStr{1},xLabelStr{2},xLabelStr{3},xLabelStr{4},...
    xLabelStr{5},xLabelStr{6},xLabelStr{7},xLabelStr{8},xLabelStr{9}});
set(h,'position', [0.03 0.03 0.43 0.015]);
set(h,'fontsize',11);
text(-173,83,'c','FontSize',7,'FontWeight','bold');
load coast;
plot(long,lat,'k-','LineWidth',0.6); hold off;
set(gca,'XLim',[-180 180],'YLim',[-90 90]);
set(gca,'XTick',-120:60:120);
set(gca,'XTickLabel',{'120W','60W','0','60E','120E'});
set(gca,'YTick',-60:30:60);
set(gca,'YTickLabel',{'60S','30S','EQ','30N','60N'});
titl = ['CTL—MODIS ',num2str2dig(ave_et_ctl_modis)];
title(titl,'fontsize',14);
grid on;

% exp-modis
subplot1(4)
v = cMin:cTick:cMax;
et_exp_modis(et_exp_modis > cMax) = cMax;
[~, h] = contourf(x,y,et_exp_modis',v);hold on;
set(h,'edgecolor','none');
axis equal;
caxis([cxisMin cxisMax]);
h = colorbar('horiz');
set(h,'XLim',[cMin cMax]);
set(h,'XTick',cMin+cTick:cTick:cMax-cTick);
set(h,'XTickLabel',{xLabelStr{1},xLabelStr{2},xLabelStr{3},xLabelStr{4},...
    xLabelStr{5},xLabelStr{6},xLabelStr{7},xLabelStr{8},xLabelStr{9}});
set(h,'position', [0.51 0.03 0.43 0.015]);
set(h,'fontsize',11);
text(-173,83,'d','FontSize',7,'FontWeight','bold');
load coast;
plot(long,lat,'k-','LineWidth',0.6); hold off;
set(gca,'XLim',[-180 180],'YLim',[-90 90]);
set(gca,'XTick',-120:60:120);
set(gca,'XTickLabel',{'120W','60W','0','60E','120E'});
set(gca,'YTick',-60:30:60);
set(gca,'YTickLabel',{'','','','',''});
titl = ['EXP—MODIS ',num2str2dig(ave_et_exp_modis)];
title(titl,'fontsize',14);
grid on;

set(gcf,'PaperPositionMode','auto')
savefigure('Revised_Figure2_ET_ctl_exp_modis.tif','PrintOptions','-dtiff', 'Size',[6.8 4.5],'LineWidth',1,'MinFontSize',8)

sysu-lipan6

不是很熟练，只想发上个代码的图

sysu-lipan6

893116547 发表于 2015-1-7 08:56
这么严重啊，那你还传，

显得悲壮一点，被发现也没事的，主要忘了删掉老师信息了

akun473

虽然不是很懂，但是赞一个

故乡月(环)

好东西，谢谢！学习一下！

1079652110

好东西