我利用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