湖北省掩膜绘制日降水填色图代码分享

LemonTree007

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import shapefile
from scipy.interpolate import griddata
from matplotlib.colors import ListedColormap, BoundaryNorm
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib.path import Path

# 设置中文字体支持
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
grid_precip = np.zeros((90,160))
# 1. 读取和处理数据
df = pd.read_csv('E:/2026//实况/20200705.csv', encoding='gbk')
df['资料时间'] = pd.to_datetime(df['资料时间'], format='%Y/%m/%d %H:%M:%S')
# 09-07
for dt in ["01", "02", "03", "04", "05", "06", "07", "08","09", "10", "11", "12","13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23"]:
    df_target = df[df['资料时间'] == pd.Timestamp('2020-07-05 '+dt+':00:00')].copy()

    # 处理异常值
    lon, lat = df_target['经度'].values, df_target['纬度'].values
    precip = df_target['过去1小时降水量'].values.copy()
    precip[precip == 999999] = np.nan
    valid_mask = ~np.isnan(precip)
    lon, lat, precip = lon[valid_mask], lat[valid_mask], precip[valid_mask]

    # 2. 插值
    grid_lon = np.arange(108.5, 116.5, 0.05)
    grid_lat = np.arange(29.0, 33.5, 0.05)
    grid_lon_mesh, grid_lat_mesh = np.meshgrid(grid_lon, grid_lat)
    grid_precip += griddata((lon, lat), precip, (grid_lon_mesh, grid_lat_mesh), method='linear', fill_value=np.nan)

# 08
df = pd.read_csv('E:/2026/LPT/实况/20200706.csv', encoding='gbk')
df['资料时间'] = pd.to_datetime(df['资料时间'], format='%Y/%m/%d %H:%M:%S')
dt = "00"
df_target = df[df['资料时间'] == pd.Timestamp('2020-07-06 ' + dt + ':00:00')].copy()

# 处理异常值
lon, lat = df_target['经度'].values, df_target['纬度'].values
precip = df_target['过去1小时降水量'].values.copy()
precip[precip == 999999] = np.nan
valid_mask = ~np.isnan(precip)
lon, lat, precip = lon[valid_mask], lat[valid_mask], precip[valid_mask]

# 2. 插值
grid_lon = np.arange(108.5, 116.5, 0.05)
grid_lat = np.arange(29.0, 33.5, 0.05)
grid_lon_mesh, grid_lat_mesh = np.meshgrid(grid_lon, grid_lat)
grid_precip += griddata((lon, lat), precip, (grid_lon_mesh, grid_lat_mesh), method='linear', fill_value=np.nan)
# print(np.max(grid_precip))
# 3. 创建掩膜
sf = shapefile.Reader("E:\\2026\地图\湖北省")
polygons = []
for shape in sf.shapes():
    for i in range(len(shape.parts)):
        start, end = shape.parts, shape.parts[i + 1] if i < len(shape.parts) - 1 else len(shape.points)
        points = shape.points[start:end]
        if len(points) > 2:
            if points[0] != points[-1]:
                points = list(points) + [points[0]]
            polygons.append(Path(points))

# 创建掩膜数组
mask = np.zeros_like(grid_lon_mesh, dtype=bool)
for i in range(grid_lon_mesh.shape[0]):
    for j in range(grid_lon_mesh.shape[1]):
        point = (grid_lon_mesh[i, j], grid_lat_mesh[i, j])
        inside = any(poly.contains_point(point) for poly in polygons)
        mask[i, j] = inside

# 应用掩膜
grid_precip_masked = grid_precip.copy()
grid_precip_masked[~mask] = np.nan

# 4. 绘图
map_crs = ccrs.PlateCarree()
fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(1, 1, 1, projection=map_crs)
ax.set_extent([108, 116.5, 29.0, 33.5])

# 绘制降水
levels = [0, 10, 25, 50, 100, 250,400,450]
colors = ['#FFFFFF', '#98FB98', '#228B22', '#1E90FF', '#0000FF','#FF00FF', '#800000']
cmap, norm = ListedColormap(colors), BoundaryNorm(levels, len(colors))
cf = ax.contourf(grid_lon_mesh, grid_lat_mesh, grid_precip_masked, levels=levels,
                     cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), alpha=0.8)

# 添加边界
for shape in sf.shapes():
    for i in range(len(shape.parts)):
        start, end = shape.parts, shape.parts[i+1] if i < len(shape.parts)-1 else len(shape.points)
        x, y = zip(*shape.points[start:end])
        if x and x[0] != x[-1]:
            x, y = list(x) + [x[0]], list(y) + [y[0]]
        ax.plot(x, y, 'k-', linewidth=1.2, transform=ccrs.PlateCarree())

ax.set_extent([108, 116.5, 29.0, 33.5], crs=ccrs.PlateCarree())
lon=[108,110,112,114,116]
lat=[29,30,31,32,33]
#坐标轴设置
ax.set_xticks(lon, crs=ccrs.PlateCarree())
ax.set_yticks(lat, crs=ccrs.PlateCarree())
ax.tick_params(labelsize=15)
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

plt.xticks(lon,fontsize=20, fontname='Times New Roman')
plt.yticks(lat,fontsize=20, fontname='Times New Roman')
plt.tick_params(axis='both', which='major', length=8,width=2)#设置坐标轴刻度的长短和粗细
# 添加颜色条
cbar = plt.colorbar(cf, orientation='horizontal', pad=0.05, aspect=50, shrink=0.7)
# 自定义刻度标签，隐藏头尾的数值
cbar.set_ticks(levels)
# 设置颜色条刻度字体
cbar.ax.tick_params(labelsize=18)
for label in cbar.ax.get_xticklabels():
    label.set_fontfamily('Times New Roman')

# 直接使用text方法，精确控制位置
cbar.ax.text(1.02, 0.4, 'mm', transform=cbar.ax.transAxes,  # y值调整到0.62
             va='center', ha='left', fontsize=18, fontfamily='Times New Roman')

plt.tight_layout()
plt.savefig('E:/2026/TP/20200705.png', dpi=300)
# plt.show()
plt.close()