我通过以下python代码,分析NC文件,并绘制了停留时间的后向轨迹图。我还需要根据山东省16地市的SHP文件(16个单独的shp文件)的面积区域,针对每个释放点,分别计算停留时间在山东省16地市区域内的在整个NC文件区域的比重,将这个比重数据计算出来,并存储在csv文件中,请帮我编写这段代码
原代码如下:
def main_analysis(nc_file, output_dir):
# ===================== 参数配置 =====================
# nc_file = "grid_time_20250921000000.nc" # 替换为您的NC文件路径
# output_dir = "./lpdm_plots2/" # 输出图片目录
height_layer = 0 # 选择高度层(0为地面层,共4层)
dpi = 300 # 图片分辨率
font()
# ===================== 数据读取 =====================
ds = nc.Dataset(nc_file)
# 提取基础变量
lons = ds.variables['longitude'][:] # (200,)
lats = ds.variables['latitude'][:] # (200,)
lon_grid, lat_grid = np.meshgrid(lons, lats) # 经纬度网格 (200x200)
lons_min = lon_grid.min()
lons_max = lon_grid.max()
lats_min = lat_grid.min()
lats_max = lat_grid.max()
# extent = [95, 130, 22, 46] # 经度范围95-130E,纬度范围22-46N
extent = [lons_min, lons_max, lats_min, lats_max]
# 提取释放点名称 (RELCOM)
relcom = ds.variables['RELCOM'][:] # shape: (14, 45)
rel_names = [relcom[i].tobytes().decode('utf-8').strip().replace('city', '') for i in range(14)]
# ------------------- 关键修改:从time变量获取起始时间 -------------------
# 读取time变量(假设为一维时间数组,单位为"hours since ..."或类似)
time_var = ds.variables['time']
if time_var is None:
raise ValueError("NC文件中未找到'time'变量,请检查文件结构")
# 获取第一个时间点的数值(通常是自基准时间的偏移量)
start_time_num = time_var[0]
end_time_num = time_var[-1]
# 将数值时间转换为datetime对象(关键步骤)
# 注意:需使用time变量的units属性(如"hours since 2025-09-21 00:00:00")
start_time = nc.num2date(
start_time_num,
units=time_var.units,
calendar=getattr(time_var, 'calendar', 'gregorian') # 默认使用公历
)
end_time = nc.num2date(
end_time_num,
units=time_var.units,
calendar=getattr(time_var, 'calendar', 'gregorian') # 默认使用公历
)
total_days = (start_time - end_time).days + 1
# -----------------------------------------------------------------------
# 提取停留时间数据 (1 ageclass, 14 pointspec, 216 time, 4 height, 200 lat, 200 lon)
spec_data = ds.variables['spec001_mr'][:]
# 提取释放点信息
rel_lons = (ds.variables['RELLNG1'][:] + ds.variables['RELLNG2'][:]) / 2
rel_lats = (ds.variables['RELLAT1'][:] + ds.variables['RELLAT2'][:]) / 2
# ===================== 绘图设置 =====================
# 地图投影设置(参考图片范围)
proj = ccrs.PlateCarree()
bounds = [
0, 0.00001,0.0005, 0.001, 0.005, 0.01, 0.1, 0.3,
0.5, 0.8, 1, 2, 3 # 10 是上界,处理 >3 的数据
]
# 2. 定义对应颜色(手动调色,严格匹配图片中的颜色分布)
colors = [
'#FFFFFF',
'#8F8FFF', # 0.0005-0.001 浅紫
'#3D3DFF', # 0.001-0.005 稍深紫
'#010FEB', # 0.005-0.01 蓝
'#03499D', # 0.01-0.1 浅蓝(接近图片深蓝)
'#06844F', # 0.1-0.3 绿
'#66D009', # 0.5-0.8 黄
'#B5E805', # 0.8-1 金黄
'#FFFF00', # 1-2 橙
'#FFA500', # 2-3 深橙
'#FF0000', # 3-10 红(覆盖 >3 区间)
'#9C0202'
]
# 3. 创建自定义配色方案和分级规范
cmap = ListedColormap(colors)
norm = BoundaryNorm(bounds, cmap.N)
# ===================== 绘制每个释放点和每个天 =====================
num_pointspec = spec_data.shape[1] # 14个释放点
for point_idx in range(num_pointspec):
logger.info(f"正在绘制{pinyin_to_hanzi[rel_names[point_idx]]}的图片")
# 此种方式为按天累加绘制
start_idx = 0
for day in range(total_days): # 9天(24h/day)
logger.info(f"{pinyin_to_hanzi[rel_names[point_idx]]}第{day+1}张图片绘制中...")
# 计算时间切片
# start_idx = day * 24 # 此种方式为按天分别绘制
end_idx = (day + 1) * 24
# 提取当前点位和当前天的停留时间数据
# shape: (24 time, 4 height, 200 lat, 200 lon)
daily_data = spec_data[0, point_idx, start_idx:end_idx, :, :, :]
# 固定高度层并积分时间维度
daily_cumulative = daily_data[:, height_layer, :, :].sum(axis=0) / 3600.0 # 转换为小时
# ------------------- 新增:生成透明度掩码(alpha数组) -------------------
# 设定透明度阈值(小于0.0005时完全透明)
alpha_threshold = 0.00001
# 生成与daily_cumulative形状相同的alpha数组(0=透明,1=不透明)
alpha = np.where(daily_cumulative < alpha_threshold, 0, 1)
# -----------------------------------------------------------------------
# 创建地图画布
fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=proj)
ax.set_extent(extent, crs=proj)
# 添加地理要素
# ------------------ 原地图地理要素 ------------------
ax.add_feature(cfeature.COASTLINE, linewidth=0.8) # 海岸线
ax.add_feature(cfeature.BORDERS, linestyle=":", linewidth=0.5) # 国界
ax.add_feature(cfeature.OCEAN, facecolor="lightblue") # 海洋
zorderidx = 1
plot_china_province(ds, ax, 'lon', 'lat', 'fig', zorderidx)
zorderidx += 1
plot_china_citys_withoutname('', ax, '', '', fig, zorderidx)
zorderidx += 1
plot_gridline('', ax, 'lon', 'lat', fig, zorderidx)
zorderidx += 1
# 绘制停留时间分布
im = ax.pcolormesh(lon_grid, lat_grid, daily_cumulative,
cmap=cmap, norm=norm, transform=proj, alpha=alpha)
# 添加颜色条
ax_height = ax.get_position().height
cbar = plt.colorbar(im, ax=ax, aspect=40, # 控制长宽比
shrink=ax_height * 0.98, # 长度
fraction=0.03,
pad=0.04,
ticks=bounds)
cbar.set_label('停留时间(小时)', size=10, labelpad=1)
cbar.ax.set_yticklabels([f'{x:.5g}' for x in bounds])
# 标记释放点位置
# ax.plot(rel_lons[point_idx], rel_lats[point_idx], 'ro', markersize=8,
# transform=proj, label='Release Point')
# 创建自定义五角星路径
star_vertices = [
(0, 1), (-0.2245, 0.3090), (-0.9511, 0.3090),
(-0.3633, -0.1180), (-0.5878, -0.8090),
(0, -0.3819), (0.5878, -0.8090),
(0.3633, -0.1180), (0.9511, 0.3090),
(0.2245, 0.3090), (0, 1)
]
star_codes = [PlotPath.MOVETO] + [PlotPath.LINETO] * 9 + [PlotPath.CLOSEPOLY]
star_path = PlotPath(star_vertices, star_codes)
# 使用自定义路径
plt.scatter(
rel_lons[point_idx], rel_lats[point_idx],
color='none', s=100, # 需要更大尺寸
marker=star_path,
edgecolor='black',
zorder=5,
transform=ccrs.PlateCarree()
)
# 调整地图边界宽度
for spine in ax.spines.values():
spine.set_linewidth(1.0)
# 确保x轴和y轴单位长度一致
ax.set_aspect('equal')
# -------------------- 新增:顶端样式文本 --------------------
# 左上角日期(示例:2025092016)
target_time = start_time - timedelta(days=day)
target_time_str = target_time.strftime("%Y%m%d%H")
path_target_time_str = target_time.strftime("%Y-%m-%d")
ax.text(0.02, 1.04, f'{target_time_str}', transform=ax.transAxes,
fontsize=14, color='black', va='top', ha='left')
# 中间地点(示例:hanzhong)
ax.text(0.5, 1.04, f'{pinyin_to_hanzi[rel_names[point_idx]]}', transform=ax.transAxes,
fontsize=14, color='black', va='top', ha='center')
# 右上角模型信息(示例:4day-LPDM)
ax.text(0.98, 1.04, f'{day + 1:02d}DAY-LPDM', transform=ax.transAxes,
fontsize=14, color='black', va='top', ha='right')
savepath = output_dir / f'{path_target_time_str}' / f'{pinyin_to_code[rel_names[point_idx]]}'
if not os.path.exists(savepath):
os.makedirs(savepath)
# 保存图片
plt.savefig(savepath / f"lpdm_{pinyin_to_code[rel_names[point_idx]]}_{path_target_time_str}.png",
dpi=dpi, bbox_inches='tight')
plt.close(fig)
logger.info(f"完成!共生成 {num_pointspec * 9} 张图片,保存至 {output_dir} 目录")
ds.close()
最新发布