python 利用cartopy绘制全国舒适指数分布图
效果图
需要用到的库
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy
import numpy as np
import netCDF4 as nc
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import matplotlib as mpl
import maskout
其中maskout是一个单独的py文件,不是用pip或anaconda安装的
绘图步骤
设置中文字体
# 设置画图字体的大小
plt.rcParams.update({'font.size': 20})
# 解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
# 解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
读取数据
file = '2022_1_1_T_dewT_u_v.nc'
nc_data_object = nc.Dataset(file)
lons = nc_data_object.variables['longitude'][:]
lats = nc_data_object.variables['latitude'][:]
t2m = (nc_data_object.variables['t2m'][:] - 273)
td = (nc_data_object.variables['d2m'][:] - 273)
w_v = numpy.sqrt(nc_data_object.variables['v10'][:] ** 2 + nc_data_object.variables['u10'][:] ** 2)
处理并计算数据
计算相对湿度
def RH1(t, td):
E0 = 6.11
E1 = E0 * 10 ** (8.6 * t / (273.3 + t))
E2 = E0 * 10 ** (8.6 * t / (273.3 + td))
RH = E2 / E1
return RH
t2m = np.sum(t2m, axis=0) / 24
td = np.sum(td, axis=0) / 24
w_v = w_v.sum(axis=0) / 24
RH = RH1(t2m, td)
I = (1.8 * t2m + 32) - 0.55 * (1 - RH) * (1.8 * t2m - 26) - 3.2 * np.sqrt(w_v)
ERA5上逐小时下载的数据,需要取平均。I为舒适度
读取地形文件并创建画布
利用cartopy.io.shapereader读取shp文件
china = shpreader.Reader('provinces.shp').geometries()
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(6, 4.5))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
设置色标
cmap = mpl.cm.RdBu_r
设置网格点并进行填图
lons_grid, lats_grid = np.meshgrid(lons, lats)
cf = ax.contourf(lons_grid, lats_grid, I, levels=50, cmap=cmap,
transform=proj, extend='both')
设置colorbar
cb = fig.colorbar(cf, shrink=0.70, orientation='horizontal', pad=0.08)
cb.set_ticks([-30, -25, -12.5, 0, 12.5, 25, 37.5, 50, 60, 80])
cb.set_ticklabels(["冷", -25, -12.5, 0, 12.5, 25, 37.5, 50, 60, "暖"])
cb.ax.set_xlabel('舒适指度', fontweight='bold')
cb.ax.tick_params(which='major', direction='in', length=3)
绘制经纬度网格线 并将上端和右侧的经纬度标注去掉
gl = ax.gridlines(alpha=0.5, linestyle='--', draw_labels=True,
dms=True, x_inline=False, y_inline=False)
gl.right_labels = 0
gl.top_labels = 0
添加海岸线(如果只要中国地图可以去掉)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=1)
绘制中国南海区域
ax_sub = fig.add_axes([0.6615, 0.3, 0.12, 0.2], projection=proj) # [*left*, *bottom*, *width*,*height*]
ax_sub.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
ax_sub.add_feature(cfeature.COASTLINE.with_scale('50m'))
china2 = shpreader.Reader('provinces.shp').geometries()
ax_sub.add_geometries(china2, ccrs.PlateCarree(), facecolor='none', edgecolor='r', linewidth=1, zorder=1)
白化(不是中国的地区不进行填色)
clip = maskout.shp2clip(cf, ax, 'china0.shp') # 白化
展示
plt.title('2022.1.1 全国舒适指度')
plt.show()
全代码
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 15 14:56:59 2021
@author: lenovo
"""
def RH1(t, td):
E0 = 6.11
E1 = E0 * 10 ** (8.6 * t / (273.3 + t))
E2 = E0 * 10 ** (8.6 * t / (273.3 + td))
RH = E2 / E1
return RH
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy
import numpy as np
import netCDF4 as nc
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import matplotlib as mpl
import maskout
# 设置画图字体的大小
plt.rcParams.update({'font.size': 20})
# 解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['SimHei']
# 解决负号乱码问题
plt.rcParams['axes.unicode_minus'] = False
file = '2022_1_1_T_dewT_u_v.nc'
nc_data_object = nc.Dataset(file)
lons = nc_data_object.variables['longitude'][:]
lats = nc_data_object.variables['latitude'][:]
t2m = (nc_data_object.variables['t2m'][:] - 273)
td = (nc_data_object.variables['d2m'][:] - 273)
w_v = numpy.sqrt(nc_data_object.variables['v10'][:] ** 2 + nc_data_object.variables['u10'][:] ** 2)
t2m = np.sum(t2m, axis=0) / 24
td = np.sum(td, axis=0) / 24
w_v = w_v.sum(axis=0) / 24
RH = RH1(t2m, td)
I = (1.8 * t2m + 32) - 0.55 * (1 - RH) * (1.8 * t2m - 26) - 3.2 * np.sqrt(w_v)
china = shpreader.Reader('provinces.shp').geometries()
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(6, 4.5))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_geometries(china, ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1, zorder=1)
cmap = mpl.cm.RdBu_r
lons_grid, lats_grid = np.meshgrid(lons, lats)
cf = ax.contourf(lons_grid, lats_grid, I, levels=50, cmap=cmap,
transform=proj, extend='both')
cb = fig.colorbar(cf, shrink=0.70, orientation='horizontal', pad=0.08)
cb.set_ticks([-30, -25, -12.5, 0, 12.5, 25, 37.5, 50, 60, 80])
cb.set_ticklabels(["冷", -25, -12.5, 0, 12.5, 25, 37.5, 50, 60, "暖"])
plt.title('2022.1.1 全国舒适指度')
cb.ax.set_xlabel('舒适指度', fontweight='bold')
cb.ax.tick_params(which='major', direction='in', length=3)
gl = ax.gridlines(alpha=0.5, linestyle='--', draw_labels=True,
dms=True, x_inline=False, y_inline=False)
gl.right_labels = 0
gl.top_labels = 0
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=1)
ax_sub = fig.add_axes([0.6615, 0.3, 0.12, 0.2], projection=proj) # [*left*, *bottom*, *width*,*height*]
ax_sub.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
ax_sub.add_feature(cfeature.COASTLINE.with_scale('50m'))
china2 = shpreader.Reader('provinces.shp').geometries()
ax_sub.add_geometries(china2, ccrs.PlateCarree(), facecolor='none', edgecolor='r', linewidth=1, zorder=1)
clip = maskout.shp2clip(cf, ax, 'china0.shp') # 白化
plt.show()