import requestsimport matplotlib.pyplot as pltimport matplotlib.dates as mdatesfrom datetime import datetimeimport numpy as np# 设置中文字体，避免图表乱码plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans', 'Arial Unicode MS']plt.rcParams['axes.unicode_minus'] = Falsedef weather_code_to_text(code): """将WMO天气代码转换为中文描述""" weather_map = { 0: "晴", 1: "晴间多云", 2: "多云", 3: "阴", 45: "雾", 48: "雾凇/冻结雾", 51: "毛毛雨：小", 53: "毛毛雨：中", 55: "毛毛雨：大", 56: "冻毛毛雨：小", 57: "冻毛毛雨：大", 61: "雨：小", 63: "雨：中", 65: "雨：大", 66: "冻雨：小", 67: "冻雨：大", 71: "雪：小", 73: "雪：中", 75: "雪：大", 77: "雪粒", 80: "阵雨：小", 81: "阵雨：中", 82: "阵雨：大/暴雨", 85: "阵雪：小", 86: "阵雪：大", 95: "雷暴：小或中", 96: "雷暴伴小冰雹", 99: "雷暴伴大冰雹", } return weather_map.get(code, f"未知({code})")def get_ecmwf_ifs_hourly(lat, lon, forecast_days=3): """获取 ECMWF IFS 模型的逐小时天气预报数据（9km 高分辨率）""" url = "https://api.open-meteo.com/v1/ecmwf" params = { "latitude": lat, "longitude": lon, "hourly": [ "temperature_2m", "precipitation", "weather_code", "wind_speed_10m", "wind_direction_10m" ], "timezone": "Asia/Shanghai", "forecast_days": forecast_days, "models": "ecmwf_ifs" } try: response = requests.get(url, params=params) response.raise_for_status() data = response.json() if "hourly" not in data: print(f"IFS 响应中未找到 'hourly' 字段，原始返回：{data}") return None hourly = data["hourly"] times = hourly.get("time", []) if not times: print("IFS 没有时间数据") return None time_objs = [] for t in times: try: time_objs.append(datetime.fromisoformat(t)) except: time_objs.append(None) result = { "time": time_objs, "temperature": hourly.get("temperature_2m", []), "precipitation": hourly.get("precipitation", []), "wind_speed": hourly.get("wind_speed_10m", []), "weather_code": hourly.get("weather_code", []), "wind_direction": hourly.get("wind_direction_10m", []) } return result except requests.exceptions.RequestException as e: print(f"IFS 网络请求错误: {e}") if hasattr(e, 'response') and e.response is not None: print(f"响应内容: {e.response.text}") return None except Exception as e: print(f"IFS 解析数据时发生错误: {e}") return Nonedef get_gfs_hourly(lat, lon, forecast_days=3): """获取 GFS (Global Forecast System) 模型的逐小时天气预报数据""" url = "https://api.open-meteo.com/v1/forecast" params = { "latitude": lat, "longitude": lon, "hourly": [ "temperature_2m", "precipitation", "weather_code", "wind_speed_10m", "wind_direction_10m" ], "timezone": "Asia/Shanghai", "forecast_days": forecast_days, "models": "gfs_seamless" } try: response = requests.get(url, params=params) response.raise_for_status() data = response.json() if "hourly" not in data: print(f"GFS 响应中未找到 'hourly' 字段，原始返回：{data}") return None hourly = data["hourly"] times = hourly.get("time", []) if not times: print("GFS 没有时间数据") return None time_objs = [] for t in times: try: time_objs.append(datetime.fromisoformat(t)) except: time_objs.append(None) result = { "time": time_objs, "temperature": hourly.get("temperature_2m", []), "precipitation": hourly.get("precipitation", []), "wind_speed": hourly.get("wind_speed_10m", []), "weather_code": hourly.get("weather_code", []), "wind_direction": hourly.get("wind_direction_10m", []) } return result except requests.exceptions.RequestException as e: print(f"GFS 网络请求错误: {e}") if hasattr(e, 'response') and e.response is not None: print(f"响应内容: {e.response.text}") return None except Exception as e: print(f"GFS 解析数据时发生错误: {e}") return Nonedef align_forecast_data(ec_data, gfs_data): """对齐两个模型的时间轴，只保留两者共同的时间点""" if not ec_data or not gfs_data: return None, None ec_time_map = {t: i for i, t in enumerate(ec_data["time"]) if t is not None} gfs_time_map = {t: i for i, t in enumerate(gfs_data["time"]) if t is not None} common_times = sorted(set(ec_time_map.keys()) & set(gfs_time_map.keys())) if not common_times: print("两个模型没有共同的时间点") return None, None aligned_ec = { "time": common_times, "temperature": [], "precipitation": [], "wind_speed": [], "weather_code": [], "wind_direction": [] } aligned_gfs = { "time": common_times, "temperature": [], "precipitation": [], "wind_speed": [], "weather_code": [], "wind_direction": [] } for t in common_times: ec_idx = ec_time_map[t] gfs_idx = gfs_time_map[t] aligned_ec["temperature"].append(ec_data["temperature"][ec_idx]) aligned_ec["precipitation"].append(ec_data["precipitation"][ec_idx]) aligned_ec["wind_speed"].append(ec_data["wind_speed"][ec_idx]) aligned_ec["weather_code"].append(ec_data["weather_code"][ec_idx]) aligned_ec["wind_direction"].append(ec_data["wind_direction"][ec_idx]) aligned_gfs["temperature"].append(gfs_data["temperature"][gfs_idx]) aligned_gfs["precipitation"].append(gfs_data["precipitation"][gfs_idx]) aligned_gfs["wind_speed"].append(gfs_data["wind_speed"][gfs_idx]) aligned_gfs["weather_code"].append(gfs_data["weather_code"][gfs_idx]) aligned_gfs["wind_direction"].append(gfs_data["wind_direction"][gfs_idx]) return aligned_ec, aligned_gfsdef print_model_data(ec_data, gfs_data, lat, lon, forecast_days): """控制台输出两个模型的对比数据（前12小时示例）""" if not ec_data or not gfs_data: print("无有效数据") return ec_times = ec_data["time"] print(f"\nECMWF IFS vs GFS 逐小时预报对比") print(f"经纬度: {lat}, {lon}") print(f"预报天数: {forecast_days}") # 显示两个模型的预报起始时间（作为起报时间参考） if ec_times and len(ec_times) > 0: ec_start = ec_times[0] gfs_start = gfs_data["time"][0] if gfs_data["time"] else None print(f"IFS 预报起始时间: {ec_start.strftime('%Y-%m-%d %H:%M')} (北京时间)") if gfs_start: print(f"GFS 预报起始时间: {gfs_start.strftime('%Y-%m-%d %H:%M')} (北京时间)") total_hours = len(ec_times) print(f"预报时效: {total_hours} 小时 (未来 {forecast_days} 天)") print("=" * 130) print(f"{'时间':<18}{'变量':<12}{'IFS':<12}{'GFS':<12}{'差值 (GFS-IFS)':<15}") print("-" * 130) max_display = min(12, len(ec_times)) variables = [ ("temperature", "气温(°C)", "temperature"), ("precipitation", "降水(mm)", "precipitation"), ("wind_speed", "风速(m/s)", "wind_speed"), ("wind_direction", "风向(°)", "wind_direction") ] for i in range(max_display): dt = ec_times[i] if dt is None: continue time_str = dt.strftime("%m-%d %H:%M") for var_key, var_name, data_key in variables: ec_val = ec_data[var_key][i] if i < len(ec_data[var_key]) else None gfs_val = gfs_data[var_key][i] if i < len(gfs_data[var_key]) else None ec_str = f"{ec_val:.1f}" if isinstance(ec_val, (int, float)) else "N/A" gfs_str = f"{gfs_val:.1f}" if isinstance(gfs_val, (int, float)) else "N/A" if isinstance(ec_val, (int, float)) and isinstance(gfs_val, (int, float)): diff = gfs_val - ec_val diff_str = f"{diff:+.1f}" else: diff_str = "N/A" print(f"{time_str:<18}{var_name:<12}{ec_str:<12}{gfs_str:<12}{diff_str:<15}") # 天气现象单独处理 ec_wcode = ec_data["weather_code"][i] if i < len(ec_data["weather_code"]) else None gfs_wcode = gfs_data["weather_code"][i] if i < len(gfs_data["weather_code"]) else None ec_weather = weather_code_to_text(ec_wcode) if ec_wcode is not None else "N/A" gfs_weather = weather_code_to_text(gfs_wcode) if gfs_wcode is not None else "N/A" print(f"{time_str:<18}{'天气现象':<12}{ec_weather:<12}{gfs_weather:<12}{'':<15}") print("-" * 130) if len(ec_times) > max_display: print(f"... 共 {len(ec_times)} 小时数据，仅显示前 {max_display} 小时")def plot_model_comparison(ec_data, gfs_data, lat, lon, forecast_days): """绘制两个模型的对比图表（气温、降水、风速、风向时间序列）""" if not ec_data or not gfs_data: print("无有效数据，无法绘图") return times = ec_data["time"] valid_indices = [i for i, t in enumerate(times) if t is not None] if not valid_indices: print("没有有效时间数据") return valid_times = [times[i] for i in valid_indices] # 提取各变量数据，并将 None 转换为 np.nan def safe_convert(data_list, indices): return [data_list[i] if i < len(data_list) and data_list[i] is not None else np.nan for i in indices] ec_temps = safe_convert(ec_data["temperature"], valid_indices) gfs_temps = safe_convert(gfs_data["temperature"], valid_indices) ec_precips = safe_convert(ec_data["precipitation"], valid_indices) gfs_precips = safe_convert(gfs_data["precipitation"], valid_indices) ec_winds = safe_convert(ec_data["wind_speed"], valid_indices) gfs_winds = safe_convert(gfs_data["wind_speed"], valid_indices) ec_wind_dir = safe_convert(ec_data["wind_direction"], valid_indices) gfs_wind_dir = safe_convert(gfs_data["wind_direction"], valid_indices) # 获取预报起始时间（作为起报时间参考） if valid_times: start_time = valid_times[0] end_time = valid_times[-1] total_hours = len(valid_times) time_info = f"IFS 起始: {start_time.strftime('%Y-%m-%d %H:%M')} | GFS 起始: {valid_times[0].strftime('%Y-%m-%d %H:%M')} | 时效: {total_hours}h" else: time_info = "无有效时间" # 创建 2×2 子图布局 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) fig.suptitle(f"ECMWF IFS vs GFS 逐小时预报对比分析\n" f"经纬度: {lat}, {lon} | 预报天数: {forecast_days}\n" f"{time_info}", fontsize=12) # 1. 气温对比 ax1 = axes[0, 0] ax1.plot(valid_times, ec_temps, label='ECMWF IFS (9km)', color='#1f77b4', linewidth=1.5, marker='.', markersize=2) ax1.plot(valid_times, gfs_temps, label='GFS (25km)', color='#ff7f0e', linewidth=1.5, marker='s', markersize=2) ax1.set_ylabel('气温 (°C)') ax1.legend(loc='upper right') ax1.grid(True, linestyle='--', alpha=0.6) ax1.set_title('气温对比') ax1.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d %H:%M')) ax1.xaxis.set_major_locator(mdates.HourLocator(interval=12)) plt.setp(ax1.xaxis.get_majorticklabels(), rotation=45, ha='right') # 2. 降水对比 ax2 = axes[0, 1] ax2.plot(valid_times, ec_precips, label='ECMWF IFS', color='#1f77b4', linewidth=1.5, marker='.', markersize=2) ax2.plot(valid_times, gfs_precips, label='GFS', color='#ff7f0e', linewidth=1.5, marker='s', markersize=2) ax2.set_ylabel('降水量 (mm)') ax2.legend(loc='upper right') ax2.grid(True, linestyle='--', alpha=0.6) ax2.set_title('降水对比') ax2.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d %H:%M')) ax2.xaxis.set_major_locator(mdates.HourLocator(interval=12)) plt.setp(ax2.xaxis.get_majorticklabels(), rotation=45, ha='right') # 3. 风速对比 ax3 = axes[1, 0] ax3.plot(valid_times, ec_winds, label='ECMWF IFS', color='#1f77b4', linewidth=1.5, marker='.', markersize=2) ax3.plot(valid_times, gfs_winds, label='GFS', color='#ff7f0e', linewidth=1.5, marker='s', markersize=2) ax3.set_ylabel('风速 (m/s)') ax3.set_xlabel('时间') ax3.legend(loc='upper right') ax3.grid(True, linestyle='--', alpha=0.6) ax3.set_title('风速对比') ax3.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d %H:%M')) ax3.xaxis.set_major_locator(mdates.HourLocator(interval=12)) plt.setp(ax3.xaxis.get_majorticklabels(), rotation=45, ha='right') # 4. 风向对比 ax4 = axes[1, 1] ax4.plot(valid_times, ec_wind_dir, label='ECMWF IFS', color='#1f77b4', linewidth=1.5, marker='.', markersize=2) ax4.plot(valid_times, gfs_wind_dir, label='GFS', color='#ff7f0e', linewidth=1.5, marker='s', markersize=2) ax4.set_ylabel('风向 (度)') ax4.set_xlabel('时间') ax4.legend(loc='upper right') ax4.grid(True, linestyle='--', alpha=0.6) ax4.set_title('风向对比 (0°=北风, 90°=东风, 180°=南风, 270°=西风)') ax4.set_ylim(0, 360) ax4.set_yticks([0, 90, 180, 270, 360]) ax4.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d %H:%M')) ax4.xaxis.set_major_locator(mdates.HourLocator(interval=12)) plt.setp(ax4.xaxis.get_majorticklabels(), rotation=45, ha='right') plt.tight_layout() plt.subplots_adjust(top=0.88) plt.show()# 主函数调用if __name__ == "__main__": lat, lon = 32.75, 118.25 # 江苏省南京市附近 forecast_days = 3 print("正在获取 ECMWF IFS 逐小时预报数据...") ec_data = get_ecmwf_ifs_hourly(lat, lon, forecast_days) print("正在获取 GFS 逐小时预报数据...") gfs_data = get_gfs_hourly(lat, lon, forecast_days) if ec_data and gfs_data: ec_aligned, gfs_aligned = align_forecast_data(ec_data, gfs_data) if ec_aligned and gfs_aligned: print_model_data(ec_aligned, gfs_aligned, lat, lon, forecast_days) plot_model_comparison(ec_aligned, gfs_aligned, lat, lon, forecast_days) else: print("数据对齐失败") else: print("获取数据失败，请检查网络或参数。")
10个月宝宝每天需要喝多少奶粉？
