Plotting data on a map 在地圖上用數(shù)據(jù)作圖

[原文地址](https://matplotlib.org/basemap/users/examples.html）
Basemap是python附加的一個可以在地圖上作圖的可視化工具。
由于我的學(xué)習(xí)路徑是通過Python for Data Analysis一書， 所以都在Jupyter notebook下進(jìn)行編譯。
為了圖省事也不是在標(biāo)準(zhǔn)python下而是用了Anaconda（一個用于科學(xué)計算的Python發(fā)行版，支持Linux, Mac, Windows系統(tǒng)，提供了包管理與環(huán)境管理的功能，可以很方便地解決多版本python并存、切換以及各種第三方包安裝問題。）在此環(huán)境下使用Basemap只需要在Anaconda Prompt中鍵入：conda install basemap; condat install pyproj 就可以下載這兩個包了。
首先介紹的是basemap中的一些實例對象：

• contour() : draw contour lines.（畫輪廓線）
• contourf() : draw filled contours.（畫填充后的輪廓線）
• imshow() : draw an image.(在地圖上畫圖）
• pcolor() : draw a pseudocolor plot.(偽色圖）
• pcolormesh() : draw a pseudocolor plot (faster version for regular meshes).
• plot() : draw lines and/or markers.(在地圖上畫線繪圖）
• scatter() : draw points with markers.（在地圖上畫散點圖）
• quiver() : draw vectors.（畫向量圖，三維就是曲面圖）
• barbs() : draw wind barbs (畫風(fēng)羽圖）
• drawgreatcircle() : draw a great circle（畫大圓航線）

在地圖上畫輪廓線

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# 常規(guī)操作， 導(dǎo)入numpy, matplotlib.pyplot和basemap
map = Basemap(projection='ortho',lat_0=45,lon_0=-100,resolution='l')
# 設(shè)置地圖正射投影點為北緯50度， 西經(jīng)100度，海岸線的分辨率為低 
# draw coastlines, country boundaries, fill continents.
map.drawcoastlines(linewidth=0.25)
# 畫出海岸線（描邊）
map.drawcountries(linewidth=0.25)
# 畫出國境線（描邊）
map.fillcontinents(color='coral',lake_color='aqua')
# 填充大陸, 大陸顏色為珊瑚色， 湖泊顏色為水色
map.drawmapboundary(fill_color='aqua')
# 畫出地圖邊界，海洋區(qū)域顏色為水色
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
# 每三十度畫出經(jīng)緯度線
nlats = 73; nlons = 145; delta = 2.*np.pi/(nlons-1)
lats = (0.5*np.pi-delta*np.indices((nlats,nlons))[0,:,:])
lons = (delta*np.indices((nlats,nlons))[1,:,:])
wave = 0.75*(np.sin(2.*lats)**8*np.cos(4.*lons))
mean = 0.5*np.cos(2.*lats)*((np.sin(2.*lats))**2 + 2.)
# 在規(guī)則網(wǎng)格上填充數(shù)據(jù)
x, y = map(lons*180./np.pi, lats*180./np.pi)
#投影到球面上
cs = map.contour(x,y,wave+mean,15,linewidths=1.5)
# 以x,y為基準(zhǔn)協(xié)調(diào)畫出wave+mean的輪廓線， 輪廓線條數(shù)為15. 參數(shù)詳情見
# https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.contour.html
plt.title('contour lines over filled continent background')
plt.show()
plt.savefig('contour lines over continent.jpg')

contour lines over continent

用填充輪廓線作的降水量圖

from mpl_toolkits.basemap import Basemap, cm
# cm(colormap)庫提供一系列彩色地圖
from netCDF4 import Dataset as NetCDFFile
import numpy as np
import matplotlib.pyplot as plt
#同上文，導(dǎo)入numpy, matplotlib.pyplot，導(dǎo)入netCDF4中的Dataset處理網(wǎng)絡(luò)通用
#數(shù)據(jù)格式（net common data form）

nc = NetCDFFile('nws_precip_conus_20061222.nc')
#首先在http://water.weather.gov/precip/中下載2006年12月22日的美國本土
#（不含阿拉斯加與夏威夷）的降水量數(shù)據(jù)
#導(dǎo)入我們需要用到的dataset， 值得注意的是該網(wǎng)站17年3月后的數(shù)據(jù)格式更新， 
#通過查詢變量名發(fā)現(xiàn)數(shù)據(jù)格式與之前有很大差異
print nc.variables.keys()
# 輸出查看數(shù)據(jù)中的變量名
prcpvar = nc.variables['amountofprecip']
data = 0.01*prcpvar[:]
latcorners = nc.variables['lat'][:]
loncorners = -nc.variables['lon'][:]
lon_0 = -nc.variables['true_lon'].getValue()
lat_0 = nc.variables['true_lat'].getValue()
# 標(biāo)準(zhǔn)化降水量與提取經(jīng)緯度參數(shù)
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
#創(chuàng)建圖像對象，設(shè)置圖像大小與軸線起始位置 
# create polar stereographic Basemap instance.
m = Basemap(projection='stere',lon_0=lon_0,lat_0=90.,lat_ts=lat_0,\
            llcrnrlat=latcorners[0],urcrnrlat=latcorners[2],\
            llcrnrlon=loncorners[0],urcrnrlon=loncorners[2],\
            rsphere=6371200.,resolution='l',area_thresh=10000)
#畫立體投影圖， 設(shè)置圖形上下左右四個邊界點經(jīng)緯度參數(shù)坐標(biāo)，中心點經(jīng)緯度參
#數(shù)坐標(biāo),定義地圖投影的球面半徑(默認(rèn)值為6370997米，近似于地球的半徑），
#分辨率以及閾值 
# 注: area_thresh = 10000 意味著面積小于10000平方公里的湖泊等對象將不被作圖
m.drawcoastlines()
m.drawstates()
m.drawcountries()
# 畫海岸線，州界， 國界線
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
#以10度為間隔畫出0度到北緯90度緯線， 并且在圖像左側(cè)設(shè)置緯線標(biāo)簽
meridians = np.arange(180.,360.,10.)
#以10度為間隔畫出西經(jīng)180度到本初子午線經(jīng)線， 并且在圖像下側(cè)設(shè)置經(jīng)線標(biāo)簽
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
ny = data.shape[0]; nx = data.shape[1]
lons, lats = m.makegrid(nx, ny)
# 經(jīng)緯線空間均勻
x, y = m(lons, lats) 
clevs = [0,1,2.5,5,7.5,10,15,20,30,40,50,70,100,150,200,250,300,400,500,600,750]
cs = m.contourf(x,y,data,clevs,cmap=cm.s3pcpn)
# 添加參數(shù)表，以x,y為基準(zhǔn)畫出data的輪廓線，等輪廓線參數(shù)為clevs，填充顏色畫出填充后的輪廓線
cbar = m.colorbar(cs,location='bottom',pad="5%")
#添加色標(biāo), 每個色標(biāo)占5%（一共20個色標(biāo)）
cbar.set_label('mm')
# 添加標(biāo)簽 單位：毫米
plt.title(prcpvar.long_name+' for period ending '+prcpvar.dateofdata)
# 添加圖像名
plt.show()
plt.savefig('24hrs rainfall of 20061222 for CONUS.jpg')

12/22/2006美國本土24小時降水量圖

繪制Argo浮標(biāo)

Argo在全球海域放置了3800個浮標(biāo)用于檢測淺層海（2000米內(nèi)）海洋的溫度與鹽度。這項計劃是的我們能夠模擬檢測淺層海的溫度，鹽度與海水流動速率。此圖我們通過數(shù)據(jù)定位Argo浮標(biāo)位置并將其繪制在地圖上。

from netCDF4 import Dataset, num2date
import time, calendar, datetime, numpy
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import urllib, os
#因為此次繪圖需要用到時間數(shù)據(jù)所以專門引入了datetime, num2date等工具， urllib用于從網(wǎng)頁提取數(shù)據(jù)（from url）
# os 這個模塊提供了一種方便的使用操作系統(tǒng)函數(shù)的方法(暫未理解)
filename, headers = urllib.urlretrieve('http://coastwatch.pfeg.noaa.gov/erddap/tabledap/apdrcArgoAll.nc?longitude,latitude,time&longitude>=0&longitude<=360&latitude>=-90&latitude<=90&time>=2017-12-07&time<=2017-12-14&distinct()')
#urllib.urlretrieve從取出網(wǎng)址對應(yīng)內(nèi)容并存放在臨時位置
dset = Dataset(filename)
# print dset 
# 打印測試：查看數(shù)據(jù)
# print dset.variables.keys
# 打印測試： 查看數(shù)據(jù)中的參數(shù)
lats = dset.variables['latitude'][:]
lons = dset.variables['longitude'][:]
time = dset.variables['time']
times = time[:]
t1 = times.min(); t2 = times.max()
#將原始數(shù)據(jù)中的經(jīng)緯度與時間提取出來
date1 = num2date(t1, units=time.units)
date2 = num2date(t2, units=time.units)
date1 = date1.strftime("%x")
date2 = date2.strftime("%x")
#將原始數(shù)據(jù)中的數(shù)字格式時間轉(zhuǎn)化為標(biāo)準(zhǔn)datatime時間并將datetime轉(zhuǎn)化為月/日/年格式
dset.close()
os.remove(filename)
#刪除之前用urlretrieve臨時存取的文件
m = Basemap(projection='hammer',lon_0=180)
#這里的hammer圖我真不知道如何翻譯...
x, y = m(lons,lats)
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
#投影到圖上， 畫地圖邊界，填充大陸圖色
m.scatter(x,y,1,marker='o',color='k')
#根據(jù)坐標(biāo)點繪制Argo浮標(biāo)散點， 第三個參數(shù)1指代散點大?。梢詫?shù)字1理解為單位半徑，數(shù)字越大散點越大），散點形狀設(shè)定為圓形，顏色設(shè)置為黑色
plt.title('Locations of %s ARGO floats active between %s and %s' %\
        (len(lats),date1,date2),fontsize=12)
# 設(shè)定圖像名稱 注： 分別提取了經(jīng)緯度參數(shù)長度（浮標(biāo)數(shù)量）,起始與末尾日期插入構(gòu)建圖像名
plt.show()

12/07/2017 ~ 12/13/2017 ARGO浮標(biāo)位置分布圖