上一期,对Python绘制气象实用地图做了比较详细的介绍,尽管已经能够满足部分需求了,但是,在实际的应用需求中,可能还是别的需求,那么,今天就手把手教大家如何绘制几个省份的白化等值线contour地图。另外,也算是对上一期进行补充,谈谈一些小技巧。
最后,对于QGIS强烈安利一波,不光它是免费的,而且跨平台,也能够完美的支持Python3.7了,能够替代大部分日常使用的ArcGIS功能,用起来不算很笨重!
目标:绘制西部的几个省份,并且mask掉其它区域,地图上支持中文,绘制经纬度网格线,附带经纬度信息。 工具:Python3.6+、ArcGIS/QGIS、Shapfile、一系列相关的Python库、测试数据
利用单独省份的Shapefile文件,制作一个shp文件包含新疆、西藏、甘肃、青海、四川,ArcGIS操作很简单不做介绍,至于QGIS我之前基本无从下手,相关的中文资料也很少,还是Google了“how to make shapefile in qgis”得到了解决方案,具体可以参考:Merge more than two Shapefile in QGIS[1],该帖子已经比较详细的做了介绍。只不过需要提醒一下,在“Merge Shapefiles”窗口中选中之前一同导入的各省份shp文件之后,窗口会奇怪的置底,需要移开当前界面才会发现其隐藏之处,不是闪退哦!再选定坐标系方案,最好和原来的shp文件一致。我在文末会提供相应的地图文件!
直接加载如下脚本,然后运行就好,代码附上:
#coding=utf-8
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) # 把温度转换为℃
# [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()
出图效果:
需要注意的地方:很多人发现输出的图片是没有经纬度的坐标信息附加在网格线两端的,怎么调都还是出不来。并且若是设置为出图显示,还会发现绘制的图怎么都挪不到最顶层。这个问题怎么解决?答案是,在脚本最顶部添加两行: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有极大好处的。
代码附上:
#coding=utf-8
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) # 把温度转换为℃
# [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://pan.baidu.com/s/1cTIk6j0E0SMpZZkdIG-_tg 密码:9g31
[1]
Python绘制气象实用地图: https://mp.weixin.qq.com/s/6oE9xpFj_xXzcGqxCctzRg
[2] Merge more than two Shapefile in QGIS: https://www.igismap.com/merge-two-shapefile-qgis/