首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >Python绘制气象实用地图[Code+Data](续)

Python绘制气象实用地图[Code+Data](续)

作者头像
MeteoAI
发布2019-08-23 17:57:43
5.1K3
发布2019-08-23 17:57:43
举报
文章被收录于专栏:MeteoAIMeteoAI

上一期,对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文件一致。我在文末会提供相应的地图文件!

第二步:如何调试程序

1.使用jupyter notebook;(推荐)

直接加载如下脚本,然后运行就好,代码附上:

#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()

出图效果:

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有极大好处的。

代码附上:

#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

References

[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/

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-08-22,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 MeteoAI 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 第一步:制作底图
  • 第二步:如何调试程序
    • 1.使用jupyter notebook;(推荐)
      • 2.直接在终端使用python xxx.py运行;
      • 地图、测试数据链接
        • References
        领券
        问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档