小編給大家分享一下Python如何使用cartopy包實(shí)現(xiàn)氣象實(shí)用地圖，希望大家閱讀完這篇文章之后都有所收獲，下面讓我們一起去探討吧！
創(chuàng)新互聯(lián)公司2013年成立，是專(zhuān)業(yè)互聯(lián)網(wǎng)技術(shù)服務(wù)公司，擁有項(xiàng)目網(wǎng)站建設(shè)、成都網(wǎng)站設(shè)計(jì)網(wǎng)站策劃，項(xiàng)目實(shí)施與項(xiàng)目整合能力。我們以讓每一個(gè)夢(mèng)想脫穎而出為使命，1280元棲霞做網(wǎng)站,已為上家服務(wù),為棲霞各地企業(yè)和個(gè)人服務(wù),聯(lián)系電話(huà):18980820575
利用單獨(dú)省份的Shapefile文件，制作一個(gè)shp文件包含新疆、西藏、甘肅、青海、四川，ArcGIS操作很簡(jiǎn)單不做介紹，至于QGIS我之前基本無(wú)從下手，相關(guān)的中文資料也很少，還是Google了“how to make shapefile in qgis”得到了解決方案，具體可以參考:Merge more than two Shapefile in QGIS[1],該帖子已經(jīng)比較詳細(xì)的做了介紹。只不過(guò)需要提醒一下，在“Merge Shapefiles”窗口中選中之前一同導(dǎo)入的各省份shp文件之后，窗口會(huì)奇怪的置底，需要移開(kāi)當(dāng)前界面才會(huì)發(fā)現(xiàn)其隱藏之處，不是閃退哦！再選定坐標(biāo)系方案，最好和原來(lái)的shp文件一致。我在文末會(huì)提供相應(yīng)的地圖文件!
1.使用jupyter notebook
直接加載如下腳本，然后運(yùn)行就好，代碼附上：
重慶分公司
重慶分公司
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
import cartopy.crs as ccrs
import shapefile
import matplotlib as mpl
import xarray as xr
from matplotlib.font_manager import FontProperties
import netCDF4 as nc
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib
ZHfont = matplotlib.font_manager.FontProperties(fname='/Users/zhpfu/Documents/fonts/SimSun.ttf')
plt.rcParams.update({'font.size': 20})
fig = plt.figure(figsize=[12,18])
ax = fig.add_subplot(111)
def basemask(cs, ax, map, shpfile):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if shape_rec.record[0] >= 0:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append(map(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform = ax.transData)
for contour in cs.collections:
contour.set_clip_path(clip)
def makedegreelabel(degreelist):
labels=[str(x)+u'°E' for x in degreelist]
return labels
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把溫度轉(zhuǎn)換為℃
# [west,east,south,north]
m = Basemap(llcrnrlon=70.,
llcrnrlat=25,
urcrnrlon=110,
urcrnrlat=50,
resolution = None,
projection = 'cyl')
cbar_kwargs = {
'orientation': 'horizontal',
'label': 'Temperature(℃)',
'ticks': np.arange(-30,30+1,10),
'pad': -0.35,
'shrink': 0.9
}
# 畫(huà)圖
levels = np.arange(-30,30+1,1)
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
m.readshapefile('west','West',color='k',linewidth=1.2)
basemask(cs, ax, m, 'west')
# draw parallels and meridians.
# label parallels on right and top
# meridians on bottom and left
parallels = np.arange(25.,50.+1,5.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[2, 3],fontsize= 14) # ha= 'right'
meridians = np.arange(70.,110.+1,5.)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[2, 3],fontsize= 14)
plt.ylabel('') #Remove the defult lat / lon label
plt.xlabel('')
plt.rcParams.update({'font.size':25})
ax.set_title(u'中國(guó)西部地區(qū)部分省份',color='blue',fontsize= 25 ,fontproperties=ZHfont) # 2m Temperature
#經(jīng)度：87.68 ， 緯度：43.77
bill0 = 87.68
tip0 = 43.77
plt.scatter(bill0, tip0,marker='.',s=120 ,color ="r",zorder=2)
#經(jīng)度：103.73 ， 緯度：36.03
bill1 = 103.73
tip1 = 36.03
plt.scatter(bill1, tip1,marker='.',s=120 ,color ="r" ,zorder=2)
#經(jīng)度：101.74 ， 緯度：36.56
bill2 = 101.74
tip2 = 36.56
plt.scatter(bill2, tip2,marker='.',s=120 ,color ="r",zorder=2 )
bill3 = 104.1
tip3 = 30.65
plt.scatter(bill3, tip3,marker='.',s=120 ,color ="r",zorder=2 )
bill4 = 91.11
tip4 = 29.97
plt.scatter(bill4, tip4,marker='.',s=120 ,color ="r",zorder=2)
plt.rcParams.update({'font.size':18})
plt.text(bill0-2.0, tip0+0.3, u"烏魯木齊",fontsize= 18 ,color ="r",fontproperties=ZHfont)
plt.text(bill1-1., tip1+0.3, u"蘭州" ,fontsize= 18 ,color ="r",fontproperties=ZHfont)
plt.text(bill2-1., tip2+0.3, u"西寧" ,fontsize= 18 ,color ="r",fontproperties=ZHfont)
plt.text(bill3-1., tip3+0.3, u"成都" ,fontsize= 18 ,color ="r",fontproperties=ZHfont)
plt.text(bill4-1., tip4+0.3, u"拉薩" ,fontsize= 18 ,color ="r",fontproperties=ZHfont)
# Save & Show figure
plt.savefig("West_China_mask.png", dpi=300, bbox_inches='tight')
plt.show()出圖效果：
2.直接在終端使用python xxx.py運(yùn)行；
需要注意的地方：很多人發(fā)現(xiàn)輸出的圖片是沒(méi)有經(jīng)緯度的坐標(biāo)信息附加在網(wǎng)格線(xiàn)兩端的，怎么調(diào)都還是出不來(lái)。并且若是設(shè)置為出圖顯示，還會(huì)發(fā)現(xiàn)繪制的圖怎么都挪不到最頂層。這個(gè)問(wèn)題怎么解決？答案是，在腳本最頂部添加兩行：import matplotlib; matplotlib.use('TkAgg')。你會(huì)發(fā)現(xiàn)看圖置頂?shù)膯?wèn)題解決了，而且網(wǎng)格線(xiàn)兩端也會(huì)正常出現(xiàn)經(jīng)緯度信息。此外，建議保存的圖片和腳本名稱(chēng)一致，解決方案：
#代碼頭部 import os,sys #代碼尾部 (filename, extension) = os.path.splitext(os.path.basename(sys.argv[0])) plt.savefig(filename+".png", dpi=600, bbox_inches='tight')
當(dāng)然，jupyter notebook是不會(huì)出現(xiàn)這個(gè)問(wèn)題的。至于什么原因大家可以去自行了解一下。還是那句話(huà)，遇到錯(cuò)誤信息了，最值得信賴(lài)的還是Google大法，學(xué)會(huì)如何使用Google，絕對(duì)是對(duì)debug有極大好處的。
代碼附上：
import matplotlib
matplotlib.use('TkAgg')
import os,sys
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
import cartopy.crs as ccrs
import shapefile
import matplotlib as mpl
import xarray as xr
from matplotlib.font_manager import FontProperties
import netCDF4 as nc
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib
ZHfont = matplotlib.font_manager.FontProperties(fname='/Users/zhpfu/Documents/fonts/SimSun.ttf')
plt.rcParams.update({'font.size': 20})
fig = plt.figure(figsize=[12,18])
ax = fig.add_subplot(111)
def basemask(cs, ax, map, shpfile):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if shape_rec.record[0] >= 0:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append(map(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform = ax.transData)
for contour in cs.collections:
contour.set_clip_path(clip)
def makedegreelabel(degreelist):
labels=[str(x)+u'°E' for x in degreelist]
return labels
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把溫度轉(zhuǎn)換為℃
# [west,east,south,north]
m = Basemap(llcrnrlon=70.,
llcrnrlat=25,
urcrnrlon=110,
urcrnrlat=50,
resolution = None,
projection = 'cyl')
cbar_kwargs = {
'orientation': 'horizontal',
'label': 'Temperature(℃)',
'ticks': np.arange(-30,30+1,10),
'pad': -0.35,
'shrink': 0.9
}
# 畫(huà)圖
levels = np.arange(-30,30+1,1)
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
m.readshapefile('west','West',color='k',linewidth=1.2)
basemask(cs, ax, m, 'west')
# draw parallels and meridians.
# label parallels on right and top
# meridians on bottom and left
parallels = np.arange(25.,50.+1,5.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[2, 3],fontsize= 14) # ha= 'right'
meridians = np.arange(70.,110.+1,5.)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[2, 3],fontsize= 14)
plt.ylabel('') #Remove the defult lat / lon label
plt.xlabel('')
plt.rcParams.update({'font.size':25})
ax.set_title(u'中國(guó)西部地區(qū)部分省份',color='blue',fontsize= 25 ,fontproperties=ZHfont) # 2m Temperature
#經(jīng)度：87.68 ， 緯度：43.77
bill0 = 87.68
tip0 = 43.77
plt.scatter(bill0, tip0,marker='.',s=120 ,color ="r",zorder=2)
#經(jīng)度：103.73 ， 緯度：36.03
bill1 = 103.73
tip1 = 36.03
plt.scatter(bill1, tip1,marker='.',s=120 ,color ="r">看完了這篇文章，相信你對(duì)“Python如何使用cartopy包實(shí)現(xiàn)氣象實(shí)用地圖”有了一定的了解，如果想了解更多相關(guān)知識(shí)，歡迎關(guān)注創(chuàng)新互聯(lián)行業(yè)資訊頻道，感謝各位的閱讀！
