python 利用cartopy绘制全国舒适指数分布图

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