第一步：制作底圖

利用單獨省份的Shapefile檔案，制作一個shp檔案包含新疆、西藏、甘肅、青海、四川，ArcGIS操作很簡單不做介紹，至于QGIS我之前基本無從下手，相關的中文資料也很少，還是Google了“how to make shapefile in qgis”得到了解決方案，具體可以參考:Merge more than two Shapefile in QGIS[1],該帖子已經比較詳細的做了介紹，只不過需要提醒一下，在“Merge Shapefiles”視窗中選中之前一同匯入的各省份shp檔案之后，視窗會奇怪的置底，需要移開當前界面才會發現其隱藏之處，不是閃退哦！再選定坐標系方案，最好和原來的shp檔案一致，我在文末會提供相應的地圖檔案!

第二步：如何除錯程式

1.使用jupyter notebook

直接加載如下腳本，然后運行就好，代碼附上：

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 = https://www.cnblogs.com/hhh188764/p/(ds['t2m'][0,::-1,:] - 273.15) # 把溫度轉換為℃

# [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
}

# 畫圖
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'中國西部地區部分省份',color='blue',fontsize= 25 ,fontproperties=ZHfont) # 2m Temperature

#經度：87.68 ， 緯度：43.77
bill0 = 87.68
tip0  =  43.77
plt.scatter(bill0, tip0,marker='.',s=120 ,color ="r",zorder=2)

#經度：103.73 ， 緯度：36.03
bill1 =  103.73
tip1  = 36.03
plt.scatter(bill1, tip1,marker='.',s=120 ,color ="r" ,zorder=2)

#經度：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運行；

需要注意的地方：很多人發現輸出的圖片是沒有經緯度的坐標資訊附加在網格線兩端的，怎么調都還是出不來，并且若是設定為出圖顯示，還會發現繪制的圖怎么都挪不到最頂層，這個問題怎么解決？答案是，在腳本最頂部添加兩行：import matplotlib; matplotlib.use('TkAgg')，你會發現看圖置頂的問題解決了，而且網格線兩端也會正常出現經緯度資訊，此外，建議保存的圖片和腳本名稱一致，解決方案：

#代碼頭部
import os,sys

#代碼尾部
(filename, extension) = os.path.splitext(os.path.basename(sys.argv[0]))
plt.savefig(filename+".png", dpi=600, bbox_inches='tight')

當然，jupyter notebook是不會出現這個問題的，至于什么原因大家可以去自行了解一下，還是那句話，遇到錯誤資訊了，最值得信賴的還是Google大法，學會如何使用Google，絕對是對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 = https://www.cnblogs.com/hhh188764/p/(ds['t2m'][0,::-1,:] - 273.15) # 把溫度轉換為℃

# [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
}

# 畫圖
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'中國西部地區部分省份',color='blue',fontsize= 25 ,fontproperties=ZHfont) # 2m Temperature

#經度：87.68 ， 緯度：43.77
bill0 = 87.68
tip0  =  43.77
plt.scatter(bill0, tip0,marker='.',s=120 ,color ="r",zorder=2)

#經度：103.73 ， 緯度：36.03
bill1 =  103.73
tip1  = 36.03
plt.scatter(bill1, tip1,marker='.',s=120 ,color ="r" ,zorder=2)

#經度：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
(filename, extension) = os.path.splitext(os.path.basename(sys.argv[0]))
plt.savefig(filename+".png", dpi=600, bbox_inches='tight')
plt.show()

轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/3032.html

標籤：Python

上一篇：Python例外處理

下一篇：Python的內置資料結構

Python資料可視化：cartopy包——氣象實用地圖

第一步：制作底圖

第二步：如何除錯程式