本帖最后由 斥鷃 于 2016-10-31 19:15 編輯

多謝@我和小云 后面的補充以及版主們的編輯,，大家可以看看一樓的推薦帖，的確比我的方法簡單很多,，在此推薦,！另裂墻推薦——有關(guān)無鋸齒白化問題可以參考@Lighting 的MATLAB對地圖進(jìn)行白化

論壇有很多人問到過這個問題，解決這個問題也有助于matlab在地學(xué)中的應(yīng)用,，前天LZ花了幾個小時看help map的文字（英語要學(xué)好?。。,?！{:soso_e134:}），粗略地研究了一下,，總算找到了幾個有用的函數(shù),，現(xiàn)在把函數(shù)和做法分享給大家，以下是一張效果圖（網(wǎng)格大小的關(guān)系，還是有些地方空了多了~）：

首先是論壇里的一些資源,，其實繪制地理類的圖形,，我一般都推薦用grads、ArcGis（不是很熟悉就不提NCL,、surfer和MeteoInfo等等啦,，論壇里關(guān)于這幾個軟件的帖子也很多，不一一舉例了）~有關(guān)這兩個軟件底圖制作和繪制在論壇里面有很詳細(xì)的教程,，大家可以參考（grads中maskout方法,，ArcGis底圖制作，ArcGis站點繪圖參考）,，而有關(guān)matlab底圖方法,，我估計m_map工具箱中會有一些解決方法，但是也不是很熟悉,，我就用matlab自帶的map工具箱中的函數(shù)進(jìn)行處理了（不清楚各版本有沒有裝,，改天看到工具箱下載的時候我共享一下）。

一,、預(yù)備知識：

1.地理參考（R）：地理參考（具體內(nèi)容可在help大部分map工具箱中的函數(shù)得到）是指定位矩陣用的矩陣,，通常為3*2的矩陣，通過公式

[lon lat] = [row col 1] * R（3*2）    (具體實例見第三部分中程序框)

對矩陣的行列坐標(biāo)與經(jīng)緯度進(jìn)行轉(zhuǎn)換（也有1*3的形式,，有興趣的朋友可以自己了解一下）,，是map工具箱處理地理數(shù)據(jù)的基礎(chǔ)；

2.矢量數(shù)據(jù)：一般矢量數(shù)據(jù)分為多邊形,、線和點數(shù)據(jù),，我們這里采用的shape文件的多邊形數(shù)據(jù)是利用線數(shù)據(jù)定義邊界的，在matlab里shape文件以結(jié)構(gòu)數(shù)組的形式保存,，其中經(jīng)緯度信息保存在X與Y屬性中；

3.2-D矩陣數(shù)據(jù)：由行標(biāo)(row)和列標(biāo)(col)來表示位置的數(shù)據(jù),，需要通過坐標(biāo)系來對數(shù)據(jù)進(jìn)行對應(yīng)進(jìn)行繪圖,，常用的繪圖參考方式有meshgrid函數(shù)生成的參考系和map中的地理參考(R)；

4.matlab中NAN（空值）的運算規(guī)則：NAN值除賦值外的任何數(shù)學(xué)運算返回值為NAN,，如NAN+1=NAN,，NAN*0=NAN；

5.在map工具箱中,，注意應(yīng)用部分函數(shù)習(xí)慣是先引用緯度值,，再引用經(jīng)度值，需要注意分別此類函數(shù),。

二,、邊界線的制作：

1.讀取文件：map=shaperead(filename)：讀取shape文件，這里建議讀取多邊形文件，通用些,。shape文件大陸部分的獲取可以由清風(fēng)提供的網(wǎng)址獲取,，不過注意邊境問題哦~另外要再做氣候區(qū)或者若干地區(qū)的底圖，用ArcGis分析工具等都可以做,，可以見ArcGis相關(guān)教程,，這里順便貢獻(xiàn)一個ArcGis的中文幫助網(wǎng)址。

2.繪制你所要繪制的圖形：這個不多說,，要畫什么畫什么,，一般是繪制等值線圖，用contour或contourf畫,，有人也許會問需不需要用contourm畫,，其實隨便，兩者的區(qū)別是前者用meshgrid生成的坐標(biāo)系來繪圖,，后者通過地理參考(R)或者經(jīng)緯度數(shù)據(jù)繪圖,，對matlab來說成圖都一樣。

3.加入邊界線：

hold on

plot(map.X,map.Y,option)

hold off

就是普通的加線了,，沒什么特別的,，option是可選的參數(shù)，美化一下繪圖而已,。

三,、掩膜的制作：

1.讀取文件，同上1,；

2.生成掩膜文件：R=makerefmat('Rastersize',size(Z),'Latlim',[lat_first lat_last],'Lonlim',[lon_first lon_last])：其中Z是你要繪圖的矩陣,，Latlim和lonlim兩個參數(shù)確定你要畫全圖的經(jīng)緯度范圍。這個函數(shù)是在后面轉(zhuǎn)換數(shù)據(jù)時候用的,，它的函數(shù)形式不止這一個（見help）,，但是我選這個用，具體原因大家可以嘗試對比一下,，因為直接使用定義經(jīng)緯度起始點和經(jīng)緯度點密度的方法生成的矩陣會出現(xiàn)少邊界點的情況,，而如上方式定義地理參考會產(chǎn)生一個縮減量使size(R)=size(Z)；

注意：鑒于有朋友提醒,，2008版前的matlab中makerefmat可能沒有'RasterSize'參數(shù)的設(shè)置,，我試了一下，以下提供上面所用方式與R=makerefmat(x0,y0,dx,dy)通用方式的轉(zhuǎn)換：

轉(zhuǎn)換要點：

1.步長dx,、dy由繪圖范圍和繪圖格點總數(shù)的比值確定,；

2.起始經(jīng)緯度為繪圖范圍的左下角經(jīng)緯度值加上半步長值。

例如：

R=makerefmat('RasterSize',size(Z),'Lonlim',[97 107],'Latlim',[21 30])

等價于

R=makerefmat(97+(107-97)/size(Z,1)/2,21+(30-21)/size(Z,2)/2,(107-97)/size(Z,1),(30-21)/size(Z,2))

>> R=makerefmat('RasterSize',size(Z),'Lonlim',[97 107],'Latlim',[21 30])

R =

0    0.0989

0.0990         0

96.9505   20.9505

>> R=makerefmat(97,21,0.1,0.1)

R =

0    0.1000

0.1000         0

96.9000   20.9000復(fù)制代碼有關(guān)投影矩陣的內(nèi)容在《經(jīng)緯度和矩陣行列標(biāo)的轉(zhuǎn)換（基于MAP工具箱）及一個副高西伸脊點程序》里還有介紹,，不過關(guān)注度有點低{:soso_e154:},，在這里打個廣告吧~

MASK=vec2mtx(map.Y,map.X,Z,R,'filled')：產(chǎn)生了一個與繪圖矩陣Z同階的,，map多邊形區(qū)域外賦值為2，多邊形邊界賦值為1,，多邊形內(nèi)賦值為0的矩陣,，同樣，vec2mtx函數(shù)也有其他輸入形式,，不過會有與上面地理參考生成時同樣的問題,，具體見help；

最后,，將MASK中值>1的區(qū)域設(shè)置為nan,，邊界設(shè)置為0，掩膜制作完成,；

3.繪圖方法：繪圖時繪制矩陣為(Z+MASK),，我這里使用的是加法方法，在2中掩膜設(shè)置方法改變一下就可以設(shè)置乘法方法的掩膜,，原理同樣是基于NAN與任何數(shù)的數(shù)學(xué)運算結(jié)果為NAN,。

附上繪制上面那幅圖的源程序，其中EOF_used是降水?dāng)?shù)據(jù)EOF的第一向量,，gy_locat是站點的站點號和經(jīng)緯度信息（測試數(shù)據(jù)見下,，數(shù)據(jù)已直接提取了插值后的Z值了，跳過了代碼的第二段）：

[lon lat]=meshgrid([97:0.1:107],[21:0.1:30]);

% Z=griddata(gy_locat(:,2),gy_locat(:,3),EOF_used(:,1),lon,lat,'v4');

yunnan=shaperead('yunnan.shp');

R=makerefmat('RasterSize',size(Z),'Lonlim',[97 107],'Latlim',[21 30]);

MASK=vec2mtx(yunnan.Y,yunnan.X,Z,R,'filled');

MASK(find(MASK>1))=nan;

MASK(find(MASK==1))=0;

contourf(lon,lat,Z+MASK,30);

shading flat

colorbar

hold on

plot(yunnan.X,yunnan.Y,'-k','linewidth',3)

hold off

不過在繪圖過程中shading flat有時候會出現(xiàn)報錯信息,，這個就不是很清楚原因了,，請各路大神指教~不過轉(zhuǎn)換方法本身應(yīng)該沒什么問題。

寫在最后：大家要是有什么好的matlab掩膜或者邊界線制作的方法都跟上一貼吧,，我也學(xué)習(xí)學(xué)習(xí),，交流一下。這個是折騰好久才弄出來的,，只是map工具箱的冰山一角,，希望大家也多有探索精神，利用身邊的資源哈~相信matlab在地學(xué)繪圖中的應(yīng)用會走得更遠(yuǎn),。

應(yīng)大家的要求,，上傳了一份處理過的測試數(shù)據(jù)：