Python实战NCDC气象数据从下载到可视化的全流程解析气象数据分析正成为环境科学、农业规划和能源管理等领域的重要工具。美国国家气候数据中心(NCDC)提供的ISD-Lite数据集因其全球覆盖和长期连续性备受研究者青睐。但对于刚接触这类数据的开发者来说从原始文件到可用分析结果之间往往存在技术鸿沟。本文将用Python构建一个完整的数据处理流水线涵盖FTP自动下载、固定宽度格式解析、数据清洗转换和可视化输出全流程。1. 环境准备与数据获取在开始处理气象数据前我们需要搭建合适的工作环境。推荐使用Anaconda创建独立的Python 3.8环境这将确保依赖管理的整洁性。核心工具链包括requests处理FTP连接与文件下载pandas数据清洗与分析的主力工具openpyxlExcel文件输出支持matplotlib基础可视化功能安装这些依赖只需一行命令pip install requests pandas openpyxl matplotlibNCDC的ISD-Lite数据通过FTP服务器公开提供文件按年份组织在固定目录结构中。例如要获取2022年中国区域的数据我们需要先确定目标气象站点的ID。这个映射关系可以在配套的isd-history.csv文件中找到import pandas as pd stations pd.read_csv(https://www1.ncdc.noaa.gov/pub/data/noaa/isd-history.csv) china_stations stations[stations[CTRY] CH].sort_values(BEGIN, ascendingFalse) print(china_stations[[USAF, STATION NAME]].head())选择目标站点后构造下载URL的格式为ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/{年}/{站点ID}-99999-{年}.gz。以下代码实现自动下载和解压import requests import gzip import shutil def download_isd_data(station_id, year, save_path): url fftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/{year}/{station_id}-99999-{year}.gz with requests.get(url, streamTrue) as r: with open(f{save_path}.gz, wb) as f: for chunk in r.iter_content(chunk_size8192): f.write(chunk) with gzip.open(f{save_path}.gz, rb) as f_in: with open(save_path, wb) as f_out: shutil.copyfileobj(f_in, f_out)2. 解析ISD-Lite固定宽度格式ISD-Lite采用固定宽度格式存储数据每行包含12个气象要素。与CSV不同这种格式需要通过精确的列位置定义来提取数据。我们先定义各字段的解析规则字段序号名称起始位置长度单位缩放因子缺失值1年份14-1-2月份62-1-3日期92-1-4小时122-1-5气温146摄氏度10-99996露点温度206摄氏度10-99997海平面气压266百帕10-99998风向326度1-99999风速386m/s10-999910云量代码446-1-9999111小时降水量506mm10-9999126小时降水量566mm10-9999使用pandas的read_fwf函数可以高效解析这种格式。我们需要特别注意处理缺失值(-9999)和单位转换def parse_isd_file(file_path): colspecs [(0,4), (5,7), (8,10), (11,13), (13,19), (19,25), (25,31), (31,37), (37,43), (43,49), (49,55), (55,61)] df pd.read_fwf(file_path, colspecscolspecs, headerNone, names[year,month,day,hour,temp,dew_point, pressure,wind_dir,wind_speed,cloud_code, precip_1h,precip_6h]) # 处理缺失值和单位转换 numeric_cols df.columns[4:] df[numeric_cols] df[numeric_cols].replace(-9999, np.nan) scale_factors {temp:10, dew_point:10, pressure:10, wind_speed:10, precip_1h:10, precip_6h:10} for col, factor in scale_factors.items(): df[col] df[col] / factor # 合并日期时间列 df[datetime] pd.to_datetime(df[[year,month,day,hour]]) return df.set_index(datetime).drop(columns[year,month,day,hour])3. 数据清洗与质量控制原始气象数据常存在各种质量问题需要系统性的清洗流程。我们首先进行基础统计检查data parse_isd_file(592870-99999-2022.txt) print(data.describe())典型的数据问题包括异常值检测气温超过50°C或低于-40°C的观测值时间连续性检查是否有缺失的时间段物理一致性露点温度不应高于气温风速突变短时间内风速剧烈变化可能指示测量错误实现自动质量控制的函数示例def quality_control(df): # 标记异常值 df[temp_qc] (df[temp] 50) (df[temp] -40) df[dew_qc] df[dew_point] df[temp] # 风速突变检测 (3小时内变化10m/s) df[wind_diff] df[wind_speed].diff().abs().rolling(3).max() df[wind_qc] df[wind_diff] 10 # 生成质量报告 qc_report pd.DataFrame({ Missing_Ratio: df.isna().mean(), Pass_QC: df.filter(like_qc).mean() }) return df, qc_report对于缺失数据我们提供多种插补策略时间序列插值适合连续缺失同期均值填充使用历史同期数据多重插补考虑变量间相关性def impute_missing(df, methodlinear): if method linear: return df.interpolate() elif method seasonal: # 使用过去7天同期均值 return df.fillna(df.groupby([df.index.hour, df.index.dayofyear]).transform(mean))4. 分析与可视化输出清洗后的数据可以输出到Excel并进行可视化。我们使用pandas的ExcelWriter创建多sheet工作簿def export_to_excel(df, output_file): with pd.ExcelWriter(output_file, engineopenpyxl) as writer: # 原始数据sheet df.to_excel(writer, sheet_nameRaw Data) # 日统计sheet daily df.resample(D).agg({ temp: [mean,max,min], precip_1h: sum, wind_speed: mean }) daily.to_excel(writer, sheet_nameDaily Stats) # 月统计sheet monthly df.resample(M).agg({ temp: mean, precip_1h: sum, wind_speed: mean }) monthly.to_excel(writer, sheet_nameMonthly Stats)可视化方面matplotlib可以生成专业图表。以下是创建气温和降水双轴图的示例import matplotlib.pyplot as plt def plot_climate(df, title): fig, ax1 plt.subplots(figsize(12,6)) # 气温曲线 color tab:red ax1.set_xlabel(Date) ax1.set_ylabel(Temperature (°C), colorcolor) ax1.plot(df.index, df[temp], colorcolor, alpha0.5) ax1.tick_params(axisy, labelcolorcolor) # 降水柱状图 ax2 ax1.twinx() color tab:blue ax2.set_ylabel(Precipitation (mm), colorcolor) ax2.bar(df.index, df[precip_1h], colorcolor, alpha0.3, width0.02) ax2.tick_params(axisy, labelcolorcolor) fig.tight_layout() plt.title(title) plt.show()对于风向数据极坐标图能更好展示分布特征def plot_wind_rose(df): from math import pi # 按16方位分组 bins np.linspace(0, 360, 17) labels [N,NNE,NE,ENE,E,ESE,SE,SSE, S,SSW,SW,WSW,W,WNW,NW,NNW] wind df.dropna(subset[wind_dir,wind_speed]) wind[dir_bin] pd.cut(wind[wind_dir], binsbins, labelslabels) # 计算各方向平均风速 rose wind.groupby(dir_bin)[wind_speed].mean().reindex(labels) # 极坐标图 ax plt.subplot(111, polarTrue) theta np.linspace(0, 2*pi, 16, endpointFalse) ax.bar(theta, rose, width2*pi/16, colorblue, alpha0.5) ax.set_theta_zero_location(N) ax.set_theta_direction(-1) plt.title(Wind Rose Diagram) plt.show()5. 自动化流水线构建将上述步骤整合为完整流水线可创建可复用的数据处理类class ISDProcessor: def __init__(self, station_id, year_range): self.station_id station_id self.year_range year_range self.raw_data [] self.clean_data None def download_all(self): for year in self.year_range: try: download_isd_data(self.station_id, year, f{self.station_id}-{year}.txt) self.raw_data.append(parse_isd_file(f{self.station_id}-{year}.txt)) except Exception as e: print(fFailed to process {year}: {str(e)}) def merge_data(self): self.clean_data pd.concat(self.raw_data).sort_index() return self.clean_data def full_pipeline(self, output_prefix): self.download_all() self.merge_data() self.clean_data, _ quality_control(self.clean_data) self.clean_data impute_missing(self.clean_data) export_to_excel(self.clean_data, f{output_prefix}.xlsx) plot_climate(self.clean_data, fStation {self.station_id}) if not self.clean_data[wind_dir].isna().all(): plot_wind_rose(self.clean_data)使用示例processor ISDProcessor(592870, range(2020,2023)) processor.full_pipeline(guangzhou_climate)6. 高级分析与扩展基础流程之外我们还可以进行更深入的分析热浪检测算法def detect_heatwaves(df, temp_thresh35, duration_days3): hot_days df[temp].resample(D).max() temp_thresh hot_periods hot_days.astype(int).groupby((~hot_days).cumsum()).cumcount() 1 return hot_periods[hot_periods duration_days]气候指标计算def climate_indices(df): years df.index.year.unique() results [] for year in years: yearly df[df.index.year year] indices { year: year, mean_temp: yearly[temp].mean(), total_precip: yearly[precip_1h].sum(), hot_days: (yearly[temp].resample(D).max() 30).sum(), dry_spell: max((yearly[precip_1h] 0).astype(int).groupby( (yearly[precip_1h] 0).cumsum()).cumcount()) } results.append(indices) return pd.DataFrame(results)与地理信息系统的集成def plot_on_map(station_id, zoom8): import folium station stations[stations[USAF] station_id].iloc[0] m folium.Map(location[station[LAT], station[LON]], zoom_startzoom) folium.Marker( [station[LAT], station[LON]], popupf{station[STATION NAME]} (ID:{station_id}) ).add_to(m) return m