[繪圖] Contour平滑化的問題請教一下
我利用m_map畫圖工具包,畫了一個太平洋的海表溫度衛星圖。但是我想要在圖中加一條
contour等溫線圖,像下圖:
http://i.imgur.com/of2OiPr.png

那一條黑線是29度的等溫線圖.
可是,我畫出來卻是像這樣子:
http://i.imgur.com/APtHGik.png

我知道是解析度太高的關係,才會變這樣.我們老師說要把資料平滑化,才可以畫出平滑的
contour圖,所以想請教各位網兄,會不會寫這一段平滑化的程式,
老師說說把九個像素的資料點(像九宮格)一樣,把這九格的資料平均,
一次移動一小格,做移動平均,這樣子的話,只會損失上,下,左,右的一排,
又可以達到把數字平均,就可以畫出平滑的contour了,謝謝,
但我,我跟本是matlab的門外漢,不知道怎麼去寫,請教各位了,百拜,感激不盡.
原始衛星資料是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
推
02/22 23:35, , 1F
02/22 23:35, 1F
推
02/24 15:11, , 2F
02/24 15:11, 2F
→
02/24 15:12, , 3F
02/24 15:12, 3F
推
02/25 23:46, , 4F
02/25 23:46, 4F
→
02/25 23:46, , 5F
02/25 23:46, 5F
MATLAB 近期熱門文章
PTT數位生活區 即時熱門文章