我利用m_map畫圖工具包,畫了一個太平洋的海表溫度衛星圖。但是我想要在圖中加一條 contour等溫線圖,像下圖: http://i.imgur.com/of2OiPr.png 那一條黑線是29度的等溫線圖. 可是,我畫出來卻是像這樣子: http://i.imgur.com/APtHGik.png 我知道是解析度太高的關係,才會變這樣.我們老師說要把資料平滑化,才可以畫出平滑的 contour圖,所以想請教各位網兄,會不會寫這一段平滑化的程式, 老師說說把九個像素的資料點（像九宮格）一樣，把這九格的資料平均， 一次移動一小格，做移動平均，這樣子的話，只會損失上，下，左，右的一排， 又可以達到把數字平均，就可以畫出平滑的ｃｏｎｔｏｕｒ了，謝謝, 但我，我跟本是ｍａｔｌａｂ的門外漢，不知道怎麼去寫，請教各位了，百拜,感激不盡. 原始衛星資料是MODIS NSST,網址是: http://oceandata.sci.gsfc.nasa.gov/MODIST/Mapped/Monthly/4km/NSST/ 程式碼如下： 下載:https://dl.dropboxusercontent.com/u/5505203/sst_l3_cut6.m clc clear all % Ocean colour data from http://seawifs.gsfc.nasa.gov/SEAWIFS.html % % Take a 4km weakly average dataset and plot a map for the Strait of % Georgia and outer coast. Note that most of this code is used % for reading in and subsetting the data. path='D:\Contour test\'; %貼上資料及程式所放的資料夾路徑 cd(path); %前往路徑的資料夾 %製作資料夾 mkdir 140E-180E-15S-15N %指定經緯度及路徑 l.path=[path '140E-180E-15S-15N\']; %(1) 140E-180E-15S-15N 最大範圍 l.lon=[130 210]; l.lat=[-20 20]; b = dir('*L3m**NSST*'); for i=1:length(b); % Note - This is probably not the most efficient way to read and % handle HDF data, but I don't usually do this... % First, get the attribute data PI=hdfinfo(b(i).name); % And write it into a structure pin=[]; for k=1:60, nm=PI.Attributes(k).Name;nm(nm==' ')='_'; if isstr(PI.Attributes(k).Value), pin=setfield(pin,nm,PI.Attributes(k).Value); else pin=setfield(pin,nm,double(PI.Attributes(k).Value)); end end; for path_j=1 % lon/lat of grid corners lon=[pin.Westernmost_Longitude:pin.Longitude_Step:pin.Easternmost_Longitude]; lat=[pin.Northernmost_Latitude:-pin.Latitude_Step:pin.Southernmost_Latitude]; if l.lon(2)>=180 lon(lon<0) = lon(lon<0)+360; lon = [lon(4321:8640) lon(1:4320)]; end % Get the indices needed for the area of interest [mn,ilt]=min(abs(lat-max(l(path_j).lat))); [mn,ilg]=min(abs(lon-min(l(path_j).lon))); ltlm=ceil(diff(l(path_j).lat)/pin.Latitude_Step); lglm=ceil(diff(l(path_j).lon)/pin.Longitude_Step); % load the subset of data needed for the map limits given if l.lon(2)>=180 P=hdfread(b(i).name,'l3m_data','Index',{[],[],[]}); P=[P(:,4321:8640) P(:,1:4320)]; P=P(ilt:(ilt+ltlm),ilg:(ilg+lglm)); else P=hdfread(b(i).name,'l3m_data','Index',{[ilt ilg],[],[ltlm+1 lglm+1]}); %+1為避免周圍出現空白 end % Convert data into log(Chla) using the equations given. Blank no-data. P=double(P); P(P==65535)=NaN; P=(pin.Slope*P+pin.Intercept); LT=lat(ilt+[0:ltlm]);LG=lon(ilg+[0:lglm]); [Plg,Plt]=meshgrid(LG,LT); % Draw the map... cd(l(path_j).path) clf; %開啟畫圖程式 m_proj('miller','lon',l(path_j).lon,'lat',l(path_j).lat); %設定投影方式 m_pcolor(Plg,Plt,P); %畫出圖形 shading flat; %去除網格 hold on m_contourf(Plg,Plt,P,[29 29],'color','k'); %畫某一特定值等值線時，[ ] 內兩數字要相同 hold off m_gshhs_i('color','k'); m_grid('tickdir','out'); % m_grid('linewi',2,'tickdir','out'); % m_grid('box','fancy','tickdir','out'); set(gca,'Clim',[20,30]); % 限定上限30下限20 h=colorbar; set(get(h,'ylabel'),'String','Temperature (^{o}C)'); set(h,'ytick',[ 20:2:30 ],... 'yticklabel',[ 20:2:30 ],... 'tickdir','out','fontsize',9); title(['MODIS SST ' datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_Start_Day,26) ' -> ' ... datestr(datenum(pin.Period_End_Year,1,0)+pin.Period_End_Day,26)],... 'fontsize',14,'fontweight','bold'); saveas(gcf,... ['SST4-' datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_Start_Day,29) '--' ... datestr(datenum(pin.Period_End_Year,1,0)+pin.Period_End_Day,29) '.png'],'png'); cd(path); end disp('done.'); end % clear all --

※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 122.116.82.166