python学习记录—— 利用再分析数据绘制天气图
2024-05-17 22:56:53
1 方法数据
1.1数据
数据采用1°*1°的NCEP/NCAR的再分析数据，数据格式为grib2。
数据下载地址：FNL1*1
说明：需要用邮箱注册账号，之后按需求下载具体日期的数据，每日4个时次，间隔6hr。
数据变量介绍：可参考
1.2 方法
平台：mac os
python及库安装管理：Anaconda
使用库：pygrib,cartopy,matplotlib,numpy
编辑器：JupyterLab
pygrib库使用方法介绍
参考公众号MeteoAI
cartopy介绍
公众号云台书使
2 绘图方法
一页4图，绘制500风场及高度场、700的风场和相对湿度、850hpa的风场及高度场、海平面气压场。
代码如下：
%matplotlib inline
#库引用
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import pygrib
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
#打开文件
path='fnl_20200327_00_00.grib2'
ds=pygrib.open(path)
#数据处理
## 500hpa 数据
hgt_500 = ds.select(name='Geopotential Height', typeOfLevel='isobaricInhPa', level=500)[0]
temp_500 = ds.select(name='Temperature', typeOfLevel='isobaricInhPa', level=500)[0]
u_500 = ds.select(name='U component of wind', typeOfLevel='isobaricInhPa', level=500)[0]
v_500= ds.select(name='V component of wind', typeOfLevel='isobaricInhPa', level=500)[0]lons=hgt_500.data()[2][0,:]
lats=hgt_500.data()[1][:,0]
hgt_500=hgt_500.data()[0]*0.1 ## 单位换算为dagpm
temp_500 = temp_500.data()[0]-273.15
u_500 = u_500.data()[0]
v_500 = v_500.data()[0]
#网格化 ，生成网格点坐标矩阵
lons, lats = np.meshgrid(lons, lats)
## 700hgt_700 = ds.select(name='Geopotential Height', typeOfLevel='isobaricInhPa', level=700)[0]
temp_700 = ds.select(name='Temperature', typeOfLevel='isobaricInhPa', level=700)[0]
u_700 = ds.select(name='U component of wind', typeOfLevel='isobaricInhPa', level=700)[0]
v_700= ds.select(name='V component of wind', typeOfLevel='isobaricInhPa', level=700)[0]
rh_700=ds.select(name='Relative humidity', typeOfLevel='isobaricInhPa', level=700)[0]
hgt_700=hgt_700.data()[0]*0.1 ## 单位换算为dagpm
temp_700 = temp_700.data()[0]-273.15
u_700 = u_700.data()[0]
v_700 = v_700.data()[0]
rh_700 = rh_700.data()[0]
## 850
hgt_850 = ds.select(name='Geopotential Height', typeOfLevel='isobaricInhPa', level=850)[0]
temp_850 = ds.select(name='Temperature', typeOfLevel='isobaricInhPa', level=850)[0]
u_850 = ds.select(name='U component of wind', typeOfLevel='isobaricInhPa', level=850)[0]
v_850= ds.select(name='V component of wind', typeOfLevel='isobaricInhPa', level=850)[0]
hgt_850=hgt_850.data()[0]*0.1 ## 单位换算为dagpm
temp_850 = temp_850.data()[0]-273.15
u_850 = u_850.data()[0]
v_850 = v_850.data()[0]## 地面
slp = ds.select(name='Pressure reduced to MSL', typeOfLevel='meanSea', level=0)[0]
slp = slp.data()[0]*0.01## 换算成hpa#绘图
#一页4图#设置投影方式，地图边界
proj = ccrs.PlateCarree() #定义投影转换，后面都会用到，不必重复输入‘ccrs.PlateCarree()’
#proj=ccrs.LambertConformal(central_longitude=120.0,central_latitude=45,standard_parallels=[40])
leftlon, rightlon, lowerlat, upperlat = (70,140,15,55)
img_extent = [leftlon, rightlon, lowerlat, upperlat]#建立画布
fig=plt.figure(figsize=(12,8))#添加第一子图
ax1 = fig.add_axes([0.1, 0.55, 0.4, 0.4],projection = proj)
#contour_map(ax1,img_extent,10)#气象要素绘制##等高线
isoheight=ax1.contour(lons,lats,hgt_500,transform=proj,levels=np.arange(500,600,4),colors='k',linewidths=1.1)#等高线标注
#ax1.clabel(cs,cs.levels[::2],colors='k',fontsize=10,fmt='%.0f')
ax1.clabel(isoheight,fontsize=7,colors='k',inline_spacing=-4,fmt='%d')
#ax1.clabel(isoheight,inline=1,colors='k',fontsize=7,fmt='%d')#等温线
isotherm=ax1.contour(lons,lats,temp_500,transform=proj,levels=range(-40,40,4),colors='r',linestyles='--',linewidths=1,alpha=0.8)
ax1.clabel(isotherm,colors='r',fontsize=7,inline_spacing=-4,fmt='%.0f')#风场ax1.barbs(lons[::4,::4],lats[::4,::4],u_500[::4,::4],v_500[::4,::4],linewidth=0.4,flagcolor='k',linestyle='-',length=5,pivot='tip',barb_increments=dict(half=2,full=4,flag=20),sizes=dict(spacing=0.15,height=0.5,width=0.12),transform=proj)##定义地理坐标标签格式
xstep, ystep = 10, 10
###set_extent需要配置相应的crs，否则出来的地图范围不准确
ax1.set_extent(img_extent,crs=proj)
#### 标注坐标轴
ax1.set_xticks(np.arange(leftlon, rightlon+xstep,xstep), crs=proj)
ax1.set_yticks(np.arange(lowerlat, upperlat+ystep,ystep), crs=proj)
# zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
ax1.set_title('(a)',loc='left')#地图设置
## Add ocean, land, rivers and lakes
#ax1.add_feature(cfeature.OCEAN.with_scale('50m'))
#ax1.add_feature(cfeature.LAND.with_scale('50m'))
#ax1.add_feature(cfeature.RIVERS.with_scale('50m'))
ax1.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)
##海岸线，50m精度
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))## 本地shp文件
china = shpreader.Reader('bou2_4l.dbf').geometries()###绘制中国国界省界九段线等等
ax1.add_geometries(china, proj,facecolor='none', edgecolor='black',zorder = 1)
###添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
sub_ax1 = fig.add_axes([0.437, 0.58, 0.07, 0.1],projection = proj)
sub_ax1.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
sub_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('bou2_4l.dbf').geometries()
sub_ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)#第二子图
ax2 = fig.add_axes([0.55, 0.55, 0.4, 0.4],projection = proj)
#湿度场
rh=ax2.contourf(lons,lats,rh_700,transform=proj,levels=np.arange(50,100,10),extend='both',cmap=plt.cm.Greens)##等高线
isoheight=ax2.contour(lons,lats,hgt_700,transform=proj,levels=np.arange(200,400,4),colors='k',linewidths=1.1)#等高线标注ax2.clabel(isoheight,fontsize=7,colors='k',inline_spacing=-4,fmt='%d')#风场
ax2.barbs(lons[::4,::4],lats[::4,::4],u_700[::4,::4],v_700[::4,::4],linewidth=0.4,flagcolor='k',linestyle='-',length=5,pivot='tip',barb_increments=dict(half=2,full=4,flag=20),sizes=dict(spacing=0.15,height=0.5,width=0.12),transform=proj)##定义地理坐标标签格式
xstep, ystep = 10, 10
###set_extent需要配置相应的crs，否则出来的地图范围不准确
ax2.set_extent(img_extent,crs=proj)
#### 标注坐标轴
ax2.set_xticks(np.arange(leftlon, rightlon+xstep,xstep), crs=proj)
ax2.set_yticks(np.arange(lowerlat, upperlat+ystep,ystep), crs=proj)
# zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)
#ax1.set_title('Station',loc='center',fontsize =15)
ax2.set_title('(b)',loc='left')#地图设置
## Add ocean, land, rivers and lakes
#ax1.add_feature(cfeature.OCEAN.with_scale('50m'))
#ax1.add_feature(cfeature.LAND.with_scale('50m'))
#ax1.add_feature(cfeature.RIVERS.with_scale('50m'))
ax2.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)
##海岸线，50m精度
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))## 本地shp文件
china = shpreader.Reader('bou2_4l.dbf').geometries()
###绘制中国国界省界九段线等等
ax2.add_geometries(china, proj,facecolor='none', edgecolor='gray',zorder = 1)
###添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
sub_ax2 = fig.add_axes([0.887, 0.58, 0.07, 0.1],projection = proj)
sub_ax2.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
sub_ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('/Users/wangzelin/2020/china_basic_map/bou2_4l.dbf').geometries()
sub_ax2.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)#第三子图
ax3 = fig.add_axes([0.1, 0.1, 0.4, 0.4],projection = proj)##等高线
isoheight=ax3.contour(lons,lats,hgt_850,transform=proj,levels=np.arange(100,200,4),colors='k',linewidths=1.1)#等高线标注ax3.clabel(isoheight,fontsize=7,colors='k',inline_spacing=-4,fmt='%d')#等温线
isotherm=ax1.contour(lons,lats,temp_500,transform=proj,levels=range(-40,40,4),colors='r',linestyles='--',linewidths=1,alpha=0.8)
ax1.clabel(isotherm,colors='r',fontsize=7,inline_spacing=-4,fmt='%.0f')#风场ax3.barbs(lons[::4,::4],lats[::4,::4],u_850[::4,::4],v_850[::4,::4],linewidth=0.4,flagcolor='k',linestyle='-',length=5,pivot='tip',barb_increments=dict(half=2,full=4,flag=20),sizes=dict(spacing=0.15,height=0.5,width=0.12),transform=proj)##定义地理坐标标签格式
xstep, ystep = 10, 10
###set_extent需要配置相应的crs，否则出来的地图范围不准确
ax3.set_extent(img_extent,crs=proj)
#### 标注坐标轴
ax3.set_xticks(np.arange(leftlon, rightlon+xstep,xstep), crs=proj)
ax3.set_yticks(np.arange(lowerlat, upperlat+ystep,ystep), crs=proj)
# zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax3.xaxis.set_major_formatter(lon_formatter)
ax3.yaxis.set_major_formatter(lat_formatter)
ax3.set_title('(c)',loc='left')#地图设置
## Add ocean, land, rivers and lakes
#ax1.add_feature(cfeature.OCEAN.with_scale('50m'))
#ax1.add_feature(cfeature.LAND.with_scale('50m'))
#ax1.add_feature(cfeature.RIVERS.with_scale('50m'))
ax3.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)
##海岸线，50m精度
ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))## 本地shp文件
china = shpreader.Reader('bou2_4l.dbf').geometries()
###绘制中国国界省界九段线等等
ax3.add_geometries(china, proj,facecolor='none', edgecolor='black',zorder = 1)
###添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
sub_ax3 = fig.add_axes([0.437, 0.13, 0.07, 0.1],projection = proj)
sub_ax3.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
sub_ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('bou2_4l.dbf').geometries()
sub_ax3.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)#第四子图
ax4 = fig.add_axes([0.55, 0.1, 0.4, 0.4],projection = proj)##等高线
c_slp=ax4.contour(lons,lats,slp,transform=proj,levels=np.arange(900,1040,2.5),colors='k',linewidths=1.1)#等高线标注ax4.clabel(c_slp,fontsize=5,colors='k',inline_spacing=-10,fmt='%.1f')##定义地理坐标标签格式
xstep, ystep = 10, 10
###set_extent需要配置相应的crs，否则出来的地图范围不准确
ax4.set_extent(img_extent,crs=proj)
#### 标注坐标轴
ax4.set_xticks(np.arange(leftlon, rightlon+xstep,xstep), crs=proj)
ax4.set_yticks(np.arange(lowerlat, upperlat+ystep,ystep), crs=proj)
# zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax4.xaxis.set_major_formatter(lon_formatter)
ax4.yaxis.set_major_formatter(lat_formatter)
ax4.set_title('(d)',loc='left')
#地图设置
## Add ocean, land, rivers and lakes
#ax1.add_feature(cfeature.OCEAN.with_scale('50m'))
#ax1.add_feature(cfeature.LAND.with_scale('50m'))
#ax1.add_feature(cfeature.RIVERS.with_scale('50m'))
ax4.add_feature(cfeature.LAKES, alpha=0.5)
#ax1.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)
##海岸线，50m精度
ax4.add_feature(cfeature.COASTLINE.with_scale('50m'))## 本地shp文件
china = shpreader.Reader('bou2_4l.dbf').geometries()###绘制中国国界省界九段线等等
ax4.add_geometries(china, proj,facecolor='none', edgecolor='black',zorder = 1)
###添加南海，实际上就是新建一个子图覆盖在之前子图的右下角
sub_ax4 = fig.add_axes([0.887, 0.13, 0.07, 0.1],projection = proj)
sub_ax4.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
sub_ax4.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('bou2_4l.dbf').geometries()
sub_ax4.add_geometries(china,ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
#存储
plt.savefig('weathermap.png',format='png',dpi=500)
效果如下：
参考学习的网站
Unidata
大佬的python气象绘图教程
python学习记录—— 利用再分析数据绘制天气图相关推荐
- python学习——进阶操作：根据数据绘制省份热力地图
输入: pro_sales.csv--省份与值 数据:province--省份列,deal--值列 输出: 中国地图.html 优化--显示省份名称 使用notepad++打开中国地图.html 搜索 ...
- python气象数据可视化学习记录1——基于ERA5数据画风场和海平面气压填色叠加图
python气象数据可视化学习记录1--基于ERA5数据画风场和海平面气压填色叠加图 1. 写在前面 2. 图片效果 3. 逐步代码解析 3.1导入库 3.2 读取NC格式数据 3.3 对数据进行加工 ...
- 小样本学习记录————利用所有数据的元学习Few-shot Text Classification with Distributional Signatures
小样本学习记录----利用所有数据的元学习Few-shot Text Classification with Distributional Signatures 在计算机视觉中,低水平的模式是可以跨学 ...
- NCEP再分析数据（FNL）Python下载
NCEP再分析数据(FNL)Python下载 楼主第一次接触WRF模型,最近在摸索该如何使用WRF和相关数据如何下载,简单记录一下. 数据网址: https://rda.ucar.edu/datase ...
- Python学习记录day6-反射、常用模块
Python学习记录day6-反射.常用模块 @(学习)[python] Python学习记录day6-反射常用模块 反射 常用模块 1 sys System-specific parameters ...
- Python学习记录（一）PIL库对于图像操作方法的简单整理
Python学习记录(一)PIL库对于图像操作方法的简单整理 首先对PIL库进行一个简单的介绍:Python图像库PIL(Python Image Library)是python的第三方图像处理库,由 ...
- Python学习记录day3
2019独角兽企业重金招聘Python工程师标准>>> Python学习记录 day3 今天是银角大王武sir讲课.先回顾了上节课所学,然后讲到了面向对象思想. set set是一个 ...
- 【Python学习记录】Numpy广播机制(broadcast)
✨ 博客主页:小小马车夫的主页 ✨ 所属专栏:Python学习记录 文章目录 一.什么是Numpy广播机制 二.Numpy广播应用 三.Numpy广播规则 一.什么是Numpy广播机制 在Numpy. ...
- python学习记录——容器篇
容器 字符串 下标(索引) # 下表也称为是索引,是一个整型数字,可以是正数,也可以是负数 # 正数下标是从0开始的,表示第一个字符,-1表示最后一个字符 my_str = 'hello'h e l ...
最新文章
- malloc(0)-malloc 0 字节
- 关于DEDECMS自定义模型当中添加自定义字段后在后台添加内容后不显示解决方案...
- 续： ［转］Oracle 表空间与数据文件
- Python -- 三元表达式（三目运算符）
- Function.prototype.bind相关知识点
- 可重组合与不相邻组合
- LightBurn（激光切割排版软件）官方中文版V1.0.04 | 激光切割排版软件哪个好
- 威金VIKing病毒专杀 ，瑞星viking专杀
- 微信小程序开屏广告实现
- Android界面 Html5还是Native，说说他们的各自的优缺点。
- 互联网时代用什么来拯救你的眼睛？
- AMBA之AHB总线
- 转让测绘资质，转让天津测绘资质
- 最好用的ios数据恢复软件：PhoneRescue for Mac
- 微信小程序自定义分享按钮
- 聊聊数据分析的权重思维：找女票身材 相貌 涵养？
- java自定义注解实例
- nginx 杂文----01
- 公众号粉丝引流裂变方式有哪些？公众号裂变涨粉有哪些方式?
- Linux中查看进程的虚拟地址空间内存布局