1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
| clear clc
collocate = load('/Users/wanzi/Desktop/tsetCTH/collocate_20230600545.txt');
MODISfile = '/Users/wanzi/Desktop/tsetCTH/MYD06_L2.A2023060.0550.061.2023132225720.hdf'; MOD_CTH = double(hdfread(MODISfile,'/cloud_top_height_1km')); MOD_CTH = double(MOD_CTH); MOD_CTH(MOD_CTH == 0 | MOD_CTH == -999)=nan;
AGRIfile_CTH = '/Users/wanzi/Desktop/tsetCTH/FY4A-_AGRI--_N_DISK_1047E_L2-_CTH-_MULT_NOM_20230301054500_20230301055959_4000M_V0001.NC'; AGRI_CTH = ncread(AGRIfile_CTH,'/CTH'); AGRI_CTH = AGRI_CTH'; AGRI_CTH(AGRI_CTH > 20000 | AGRI_CTH == -999) = nan;
AGRIfile_CLP = '/Users/wanzi/Desktop/tsetCTH/FY4A-_AGRI--_N_DISK_1047E_L2-_CLP-_MULT_NOM_20230301054500_20230301055959_4000M_V0001.NC'; AGRI_CLP0 = ncread(AGRIfile_CLP,'/CLP'); AGRI_CLP0 = AGRI_CLP0'; AGRI_CLP0(AGRI_CLP0 >= 126) = nan; [m,n]=size(AGRI_CLP0); AGRI_CLP=zeros(m,n); for i=1:m for j=1:n if AGRI_CLP0(i,j) == 0 AGRI_CLP(i,j) = 0; end if AGRI_CLP0(i,j) == 1 || AGRI_CLP0(i,j) == 2 AGRI_CLP(i,j) = 1; end if AGRI_CLP0(i,j) == 3 || AGRI_CLP0(i,j) == 4 AGRI_CLP(i,j) = 2; end if AGRI_CLP0(i,j) == 5 AGRI_CLP(i,j) = 3; end end end
lines1 = 1; lines2 = length(collocate); Lat = collocate(lines1:lines2,5); Lon = collocate(lines1:lines2,6); loc = collocate(lines1:lines2,1:4); CLP_collocateMODIS = collocate(lines1:lines2,12); CTH_collocateMODIS = zeros(lines2-lines1+1,1); CLP_collocateAGRI = zeros(lines2-lines1+1,1); CTH_collocateAGRI = zeros(lines2-lines1+1,1);
for i=1:lines2-lines1+1 CLP_collocateAGRI(i,1) = AGRI_CLP(loc(i,1)+1,loc(i,2)+1); CTH_collocateAGRI(i,1) = AGRI_CTH(loc(i,1)+1,loc(i,2)+1); CTH_collocateMODIS(i,1) = MOD_CTH(loc(i,3)+1,loc(i,4)+1); end
CLP_collocateMODIS(CLP_collocateMODIS(:,1)==2) = 1; CLP_collocateMODIS(CLP_collocateMODIS(:,1)==3) = 2;
MOD_info= hdfinfo(MODISfile); MOD_COT = double(hdfread(MODISfile,'/Cloud_Optical_Thickness')); MOD_CER = double(hdfread(MODISfile,'/Cloud_Effective_Radius')); MOD_COT(MOD_COT == -9999) = 0; MOD_CER(MOD_CER == -9999) = 0; M_COT_scale = MOD_info.Vgroup.Vgroup(2).SDS(73).Attributes(5).Value; M_COT_offset= MOD_info.Vgroup.Vgroup(2).SDS(73).Attributes(6).Value; M_CER_scale = MOD_info.Vgroup.Vgroup(2).SDS(67).Attributes(5).Value; M_CER_offset= MOD_info.Vgroup.Vgroup(2).SDS(67).Attributes(6).Value; MODIS_COT = MOD_COT.*M_COT_scale+M_COT_offset; MODIS_CER = MOD_CER.*M_CER_scale+M_CER_offset;
COT_collocateMODIS = zeros(lines2-lines1+1,1); COT = collocate(lines1:lines2,20); for i=1:lines2-lines1+1 COT_collocateMODIS(i,1) = MODIS_COT(loc(i,3)+1,loc(i,4)+1); end figure scatter(COT_collocateMODIS,COT)
CLP_collocateMODIS_all = [CLP_collocateMODIS_all;CLP_collocateMODIS]; CLP_collocateAGRI_all = [CLP_collocateAGRI_all;CLP_collocateAGRI]; CTH_collocateAGRI_all = [CTH_collocateAGRI_all;CTH_collocateAGRI]; CTH_collocateMODIS_all = [CTH_collocateMODIS_all;CTH_collocateMODIS];
set(gcf,'units','pixels','position',[50,50,900,800]);
position_ax1 = [0.06,0.55,0.35,0.39]; colorbar_ax1 = [0.43,0.55,0.02,0.39]; ax1 = subplot('Position',position_ax1); scatter(Lon,Lat,7,CLP_collocateAGRI,'filled') colormap(ax1,[1 1 1;0 0 1;0.3 0.75 0.93]); caxis(ax1,[-0.5 2.5]); hcb = colorbar('position',colorbar_ax1); set(hcb,'ytick',-1:1:3,'tickdir','out','YTickLabel',{'','Clear','Water','Ice',''},'FontName','Arial','FontSize',15), set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out'); set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out'); xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ... '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ... '160°E','165°E','170°E','175°E','180°E'}) yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ... '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'}) box on title('(a) AGRI CPH') set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);
position_ax4 = [0.06,0.06,0.35,0.39]; colorbar_ax4 = [0.43,0.06,0.02,0.39]; ax4 = subplot('Position',position_ax4); scatter(Lon,Lat,7,CLP_collocateMODIS,'filled') colormap(ax4,[1 1 1;0 0 1;0.3 0.75 0.93]); caxis(ax4,[-0.5 2.5]); hcb = colorbar('position',colorbar_ax4); set(hcb,'ytick',-1:1:3,'tickdir','out','YTickLabel',{'','Clear','Water','Ice',''},'FontName','Arial','FontSize',15), set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out'); set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out'); xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ... '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ... '160°E','165°E','170°E','175°E','180°E'}) yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ... '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'}) box on title('(c) MODIS CPH') set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);
INDEX_CTH = 4; INTERVAL = 1; position_ax2 = [0.57,0.55,0.35,0.39]; colorbar_ax2 = [0.94,0.55,0.02,0.39]; ax2 = subplot('Position',position_ax2); scatter(Lon,Lat,7,CTH_collocateAGRI*0.001,'filled')
hcb=colorbar('position',colorbar_ax2); caxis(ax2,[0 INDEX_CTH]) set(hcb,'ytick',0:INTERVAL:20,'tickdir','out','FontName','Arial','FontSize',15) set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out'); set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out'); xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ... '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ... '160°E','165°E','170°E','175°E','180°E'}) yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ... '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'}) box on title('(b) AGRI CTH (km)') set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);
position_ax5 = [0.57,0.06,0.35,0.39]; colorbar_ax5 = [0.94,0.06,0.02,0.39]; ax5 = subplot('position',position_ax5); scatter(Lon,Lat,7,CTH_collocateMODIS*0.001,'filled')
hcb=colorbar('position',colorbar_ax5); caxis(ax5,[0 INDEX_CTH]) set(hcb,'ytick',0:INTERVAL:20,'tickdir','out','FontName','Arial','FontSize',15) set(gca,'YLim',[Latmin-0.5 Latmax+0.5],'YTick',-50:5:50,'tickdir','out'); set(gca,'XLim',[Lonmin-0.5 Lonmax+0.5],'XTick',40:5:180,'tickdir','out'); xticklabels({'40°E','45°E','50°E','55°E','60°E','65°E','70°E','75°E','80°E','85°E','90°E','95°E', ... '100°E','105°E','110°E','115°E','120°E','125°E','130°E','135°E','140°E','145°E','150°E','155°E', ... '160°E','165°E','170°E','175°E','180°E'}) yticklabels({'50°S','45°S','40°S','35°S','30°S','25°S','20°S','15°S','10°S','5°S','0°', ... '5°N','10°N','15°N','20°N','25°N','30°N','35°N','40°N','45°N','50°N'}) box on title('(d) MODIS CTH (km)') set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7);
clear load('./testCTH.mat') CTH_collocateMODIS_all(CLP_collocateMODIS_all ~= CLP_collocateAGRI_all,:)=[]; CTH_collocateAGRI_all(CLP_collocateMODIS_all ~= CLP_collocateAGRI_all,:)=[];
pixel_num=length(CTH_collocateMODIS_all); n=100; xmin = 0; xmax = 18; ymin = 0; ymax = 18; x = CTH_collocateMODIS_all*0.001; y = CTH_collocateAGRI_all*0.001; [ Nbin ] = Prob( x, y, xmin, xmax, ymin, ymax, n); pdata = Nbin'./pixel_num; max=max(max(pdata)); min=min(min(pdata)); pdata1 = (pdata - min)/(max-min); pdata1(pdata1==0)=nan; x_lable=xmin:0.18:xmax; y_lable=ymin:0.18:ymax;
pdata1(pdata1<0)=nan; pcolor(x_lable(1,1:end-1),y_lable(1,1:end-1),pdata1) shading flat set(gca,'FontName','Arial','FontSize',15,'LineWidth',1.7); set(gca,'xlim',[0 18],'xtick',0:3:50,'tickdir','out'); set(gca,'ylim',[0 18],'ytick',0:3:50,'tickdir','out'); xlabel('Aqua/MODIS CTH (km)') ylabel('FY-4A/AGRI CTH (km)') axis square hcb=colorbar; caxis([0 1]) set(hcb,'ytick',0:0.2:1,'tickdir','out') set(hcb,'TickLabels',{'0.0','0.2','0.4','0.6','0.8','1.0'}) hold on x1 = [0 100]; y1 = [0 100]; plot(x1,y1,'--','color','black','LineWidth',1.7);
|