专栏首页MeteoAIPython绘制气象实用地图[Code+Data](续)

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

上一期,对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/

本文分享自微信公众号 - MeteoAI(meteoai),作者:gavin7675

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-08-22

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 如何用python画图——带你入门matplotlib

    matplotlib是python中常用的一个可视化库,大多数的操作与MATLAB非常类似,所以对于从MATLAB迁移到python的朋友是非常友好的。matp...

    zhangqibot
  • Deecamp 夏令营 AI 降水预测总结

    本文作者是气科院2020届硕士生方祖亮同学,这篇文章是他参加Deecamp夏令营的一个总结。方祖亮同学本科毕业于兰州大学,目前在气科院读研三,师从俞小鼎和王秀明...

    zhangqibot
  • Python中关于底图的操作

    Python基于其强大的功能越来越成为了科学利器,气象上对精细化的要求越来越高,对于底图的制作也越来越高。本人气象出身,长期用NCL画图,但是NCL对于精细化底...

    zhangqibot
  • Java 命令行运行参数大全

    javac 用法:javac <选项> <源文件> 其中,可能的选项包括:   -g                                   ...

    汤高
  • Mac下通过VMware Fusion安装centos虚拟机操作记录

    下面介绍下利用VMware Fusion工具在Mac上安装centos虚拟机的做法: 1)下载VMware Fusion工具 下载地址(包括注册码):http:...

    洗尽了浮华
  • 数据结构学习-python实现01--0401

    经过近两年多的转行自学,乱七八糟的学了不少的东西,依然没有走到自己想要去的方向,继续学习,努力吧!

    到不了的都叫做远方
  • 为什么 HTTPS 是安全的?

    HTTP 协议是通过客户端和服务器的请求应答来进行通讯,目前协议由之前的 RFC 2616 拆分成立六个单独的协议说明(RFC 7230、RFC 7231、RF...

    民工哥
  • B2B的5种经营模式

      这是一种以买家为中心,专门为某一家公司所设计的采购型网站,它是由买方自己投资建设的。例如像英特尔、沃玛特、IBM、通用汽车、戴尔电脑等。

    用户2192970
  • python+OpenCV3 实现人脸检测 实时读取视频流

    XML文件:https://github.com/opencv/opencv/blob/master/data/haarcascades/haarcascade...

    种花家的奋斗兔
  • Android加载Assets目录中Xml布局文件

    最近由于项目开发使用到了动态布局,因为打包sdk ,sdk 这块activity 需要一些layout 文件 。而做过sdk 开发的小伙伴应该知道,layout...

    砸漏

扫码关注云+社区

领取腾讯云代金券