ARTICLE DETAIL

建站实战干货

来自一线的建站与推广经验沉淀,每一条都经过真实交付验证。

气象点面数据融合:WGS84下点栅格一致性校验与NetCDF构建

2026/10/3 8:54:14 拓冰建站 浏览量
气象点面数据融合:WGS84下点栅格一致性校验与NetCDF构建 简介本资源提供中国2020年全国范围高精度年均气温空间数据适用于GIS空间分析、气候研究、环境建模及地理信息教学等场景特别适合ArcGIS初学者与科研人员开展温度空间插值、站点-栅格协同分析或区域气候特征可视化实践。压缩包共10个文件包含核心的.shp矢量站点数据与.tif栅格影像含.prj坐标定义、.tfw地理配准、.ovr金字塔、.xml元数据及.dbf属性表等完整GIS支持文件结构规范开箱即用总大小仅100KB轻量高效。已有1637人学习下载表明其在教学演示与快速验证类项目中具备较高实用价值。用户可直接加载shp点位与tif栅格进行叠加分析完成空间匹配、统计摘要、制图输出等典型GIS操作无需额外处理即可支撑课程实验、毕业设计或科研预研中的基础气候数据分析任务。1. 中国2020年均气温数据点加栅格.zip不是“下载即用”的气象包而是空间数据融合的典型入口你双击解压这个 ZIP 包看到points.csv和grid_0.1deg.tif两个文件时第一反应可能是——“终于有现成的全国气温数据了”。但很快就会卡在CSV 里的经纬度怎么和 GeoTIFF 对不上为什么用 QGIS 打开栅格后颜色发灰、数值全为 0ArcGIS 提示“坐标系缺失”却找不到.prj文件这不是数据损坏而是中国气象观测数据落地时最常被忽略的空间基准统一问题2020 年全国 2400 国家级气象站实测均温点数据与 CMORPH 或 CMA-LSM 生成的 0.1°×0.1° 空间插值栅格面数据天然存在观测尺度差异、投影系统错位、时间代表性和单位不一致四大断层。本篇不讲气候学理论只聚焦一线工程师拿到这个 ZIP 后72 小时内完成“点面一致性校验→坐标系强制对齐→双源数据空间叠加→导出可直接喂给 PyTorch 气候模型训练的 NetCDF”这一完整链路。适合正在做城市热岛分析、农业物候建模或气象深度学习输入预处理的从业者——尤其当你发现模型在验证集上 R² 突然掉 0.3而根源只是points.csv里lon列用的是 WGS84 经度grid_0.1deg.tif却是 CGCS2000 / Albers 投影下的米制坐标。2. 解包后第一步识别数据本质拒绝“文件名即真相”这个 ZIP 的命名极具迷惑性。“中国2020年均气温”听起来像一个标准产品但实际是两类异构数据的物理打包离散观测点point与连续空间场raster。二者来源、精度、误差特征完全不同。不先拆解后续所有操作都是空中楼阁。2.1 点数据points.csv的真实结构与陷阱用pandas快速探查import pandas as pd df pd.read_csv(points.csv, encodinggbk) # 注意国产气象数据常用 GBK 编码UTF-8 会乱码 print(df.head()) print(df.info())提示若报UnicodeDecodeError99% 是编码问题。不要强行errorsignore那会 silently 破坏经纬度数值。优先试gbk、gb2312、gb18030。典型输出station_id lon lat temp_2020 province 0 50135 116.4 39.9 12.85 北京 1 50136 116.5 40.0 12.72 北京 2 50137 117.2 39.1 13.01 天津 ...关键字段解读station_id中国气象数据网CMDN标准台站编号5 位数字可反查台站元数据如海拔、建站年份lon/latWGS84 坐标系下的经纬度十进制度这是中国地面观测数据的法定发布标准但注意部分老站数据可能含 0.01° 级别录入误差temp_2020该站 2020 年日均温的算术平均值单位℃非插值结果是实测值province省级行政区名称中文字符非 ISO 代码用于后续空间聚合但需注意“内蒙古”“宁夏”等带“自治区”后缀的写法一致性2.2 栅格数据grid_0.1deg.tif的元信息深挖不要依赖文件名判断分辨率或坐标系。用rasterio直接读取元数据import rasterio with rasterio.open(grid_0.1deg.tif) as src: print(CRS:, src.crs) # 输出坐标系定义 print(Bounds:, src.bounds) # 输出地理范围左下/右上经纬度 print(Res:, src.res) # 输出像元分辨率经度方向, 纬度方向 print(Shape:, src.shape) # 输出行列数 data src.read(1) # 读取第 1 波段气温值 print(Min/Max:, data.min(), data.max())常见结果CRS: EPSG:4326 Bounds: BoundingBox(left73.0, bottom18.0, right135.0, top54.0) Res: (0.1, 0.1) Shape: (360, 620) Min/Max: -52.3 34.8⚠️ 注意EPSG:4326表示 WGS84 地理坐标系经纬度res(0.1, 0.1)表明是 0.1°×0.1° 规则网格但“0.1°”在赤道约 11km在北纬 50° 仅约 7km——这是后续空间匹配误差的物理根源。Bounds显示覆盖中国全境73°E–135°E, 18°N–54°N但需警惕栅格边缘常含 NoData 值如 -9999必须用src.nodata确认填充值。2.3 为什么必须做“点面一致性校验”点数据是 2400 个离散位置的实测值栅格是 360×620223,200 个像元的插值场。二者统计口径不同点数据每个站代表其周边 1km² 内的局地气候受地形、城市热岛、仪器高度影响显著栅格数据基于克里金插值或机器学习回归生成平滑了局部异常但可能低估山地/海岸线温度梯度不做校验就叠加等于把“体温计读数”和“红外热成像图”强行比对——数值量级可能一致但空间语义完全错位。校验目标不是让二者数值相等而是确认当points.csv中某站位于grid_0.1deg.tif的某个像元中心 5km 内时其temp_2020与该像元值的绝对偏差是否在 ±1.5℃ 内中国气象行业默认可接受插值误差阈值。这一步直接决定后续建模的物理可信度。3. 坐标系强制对齐WGS84 下的“经纬度对齐”不是万能解药很多教程说“都是 WGS84直接 overlay 就行”。这是最大误区。points.csv的lon/lat是点坐标grid_0.1deg.tif的EPSG:4326是地理坐标系但栅格的“地理坐标系”本质是经纬度网格其像元中心坐标需按球面距离计算而非平面直角坐标。当你要提取某点落入哪个像元时必须用地理空间算法而非简单四舍五入。3.1 点数据转 GeoDataFrame添加几何列并验证 CRSimport geopandas as gpd from shapely.geometry import Point # 构造 geometry 列 geometry [Point(xy) for xy in zip(df[lon], df[lat])] gdf gpd.GeoDataFrame(df, geometrygeometry, crsEPSG:4326) # 验证检查是否有坐标越界如 lon180.1 或 lat99.9 invalid gdf[~gdf.geometry.is_valid] if len(invalid) 0: print(发现无效坐标点, invalid[[station_id, lon, lat]]) # 常见修复lon 超 180→减360lat 超 90→取绝对值或丢弃 gdf gdf[gdf.geometry.is_valid] print(点数据 CRS:, gdf.crs) # 应输出 EPSG:43263.2 栅格数据重采样到统一地理网格为何必须做grid_0.1deg.tif的Bounds是(73.0, 18.0, 135.0, 54.0)但points.csv中的站点lon/lat可能落在72.99或135.01—— 这些点在原始栅格中属于 NoData 区域。直接rasterio.sample()会返回nan导致 200 个站点丢失。解决方案用rasterio.warp.reproject将栅格扩展 0.05° 边界并确保像元中心严格对齐 WGS84 十进制度网格import numpy as np from rasterio.warp import reproject, Resampling from rasterio.transform import from_origin # 定义新网格以 points.csv 的 min/max 为界扩展 0.05°分辨率保持 0.1° lon_min, lon_max gdf.geometry.x.min() - 0.05, gdf.geometry.x.max() 0.05 lat_min, lat_max gdf.geometry.y.min() - 0.05, gdf.geometry.y.max() 0.05 # 计算新像元数向上取整保证覆盖 width int(np.ceil((lon_max - lon_min) / 0.1)) height int(np.ceil((lat_max - lat_min) / 0.1)) # 构建新 transform左上角为 (lon_min, lat_max)分辨率 (0.1, -0.1) transform from_origin(lon_min, lat_max, 0.1, 0.1) # 读取原栅格 with rasterio.open(grid_0.1deg.tif) as src: # 创建新空数组 dst_array np.zeros((height, width), dtypesrc.dtypes[0]) # 重投影关键使用 bilinear 插值避免 nearest 导致阶梯效应 reproject( sourcerasterio.band(src, 1), destinationdst_array, src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crssrc.crs, # 保持 EPSG:4326 resamplingResampling.bilinear, num_threads4 ) # 保存新栅格 profile src.profile.copy() profile.update({ height: height, width: width, transform: transform, driver: GTiff, nodata: src.nodata }) with rasterio.open(grid_aligned.tif, w, **profile) as dst: dst.write(dst_array, 1)参数说明Resampling.bilinear是必须项——nearest会把 0.1° 栅格变成“马赛克”cubic过度平滑。num_threads4加速但超过 CPU 核数无益。transform中lat方向步长为负因为 GeoTIFF 的transform定义 y 轴向下为正而地理坐标系 y 轴向上为正。3.3 提取点位对应栅格值地理空间最近邻 vs 像元中心匹配错误做法int((lon - 73.0) / 0.1)直接索引——忽略地球曲率高纬度偏移可达 3km。正确做法用rasterio.sample它内部调用 GDAL 的地理空间采样器from rasterio.sample import sample_gen # 读取对齐后的栅格 with rasterio.open(grid_aligned.tif) as src: # 提取每个点的栅格值 coords [(x, y) for x, y in zip(gdf.geometry.x, gdf.geometry.y)] samples list(sample_gen(src, coords)) # 返回生成器需 list() 强制执行 # 赋值回 GeoDataFrame gdf[grid_temp] [s[0] if s[0] ! src.nodata else np.nan for s in samples] # 查看偏差统计 gdf[abs_error] (gdf[temp_2020] - gdf[grid_temp]).abs() print(偏差统计℃) print(gdf[abs_error].describe())此时你会看到大部分点误差 1.0℃但青藏高原边缘、黑龙江北部站点误差 2.5℃——这正是需要人工核查或剔除的“异常点”。4. 点面空间叠加实战生成可用于深度学习的 NetCDF 数据集模型训练需要结构化、自描述、支持多维的格式。CSV TIFF 组合无法直接喂给 PyTorch DataLoader。必须合成单一时空 NetCDF 文件包含lat,lon,temp_point,temp_grid,error四个变量。4.1 构建 NetCDF 的维度与坐标系统NetCDF 要求显式声明维度dimension和坐标变量coordinate variable。这里我们以点数据为基准构建n_points维度import netCDF4 as nc import numpy as np # 创建新 NetCDF 文件 ds nc.Dataset(china_temp_2020.nc, w, formatNETCDF4) # 定义维度 n_points len(gdf) ds.createDimension(point, n_points) ds.createDimension(time, 1) # 2020 年单一时次 # 创建坐标变量 lat_var ds.createVariable(lat, f4, (point,)) lat_var.units degrees_north lat_var[:] gdf.geometry.y.values lon_var ds.createVariable(lon, f4, (point,)) lon_var.units degrees_east lon_var[:] gdf.geometry.x.values time_var ds.createVariable(time, i4, (time,)) time_var.units days since 2020-01-01 time_var.calendar gregorian time_var[:] 0 # 2020-01-01 # 创建数据变量 temp_point ds.createVariable(temp_point, f4, (point,)) temp_point.units degree_C temp_point.long_name Observed annual mean temperature at station temp_point[:] gdf[temp_2020].values temp_grid ds.createVariable(temp_grid, f4, (point,)) temp_grid.units degree_C temp_grid.long_name Interpolated annual mean temperature from 0.1deg grid temp_grid[:] gdf[grid_temp].values error ds.createVariable(error, f4, (point,)) error.units degree_C error.long_name Absolute difference between observed and interpolated temperature error[:] gdf[abs_error].values # 添加全局属性 ds.description China 2020 annual mean temperature: station observations and 0.1deg interpolated grid values ds.history Created on str(pd.Timestamp.now()) ds.source CMA (China Meteorological Administration) CMA-LSM reanalysis ds.close()逻辑说明NetCDF 不存储几何对象所以lat/lon作为一维坐标变量temp_point等作为一维数据变量隐含“第 i 个点的 lat/lon/temp 对应同一物理位置”。这种结构被xarray.open_dataset()完美支持可直接ds.temp_point.plot.hist()可视化。4.2 添加空间权重解决“站点分布不均”导致的模型偏差中国气象站密度差异极大东部平原每万 km² 有 10 站青藏高原每万 km² 不足 0.1 站。若直接用temp_point训练模型会严重过拟合东部数据。必须引入空间权重spatial weight# 计算每个站点的 Voronoi 多边形面积泰森多边形作为权重 from scipy.spatial import Voronoi, voronoi_plot_2d import shapely.ops as ops from shapely.geometry import Polygon, MultiPolygon # 获取中国国界简化版用于裁剪 # 实际项目中应下载 GADM 或 Natural Earth 的 china_adm0.shp # 此处用 bbox 近似73°E-135°E, 18°N-54°N china_bbox Polygon([(73, 18), (135, 18), (135, 54), (73, 54)]) # 构造 Voronoi points np.column_stack([gdf.geometry.x, gdf.geometry.y]) vor Voronoi(points) # 提取 Voronoi 多边形并裁剪到中国范围 regions [] for region in vor.regions: if not -1 in region and len(region) 0: polygon Polygon([vor.vertices[i] for i in region]) clipped polygon.intersection(china_bbox) if not clipped.is_empty: regions.append(clipped.area) # 赋权面积越大权重越小稀疏区站点更珍贵 weights 1.0 / np.array(regions) weights weights / weights.sum() # 归一化 gdf[spatial_weight] weights # 写入 NetCDF weight_var ds.createVariable(spatial_weight, f4, (point,)) weight_var.long_name Spatial weight based on Voronoi polygon area weight_var[:] weights参数说明Voronoi在边界会产生无限区域含-1索引必须过滤。china_bbox是粗略矩形实际应用建议用gpd.read_file(china_boundary.shp)精确裁剪。权重归一化确保sum(weights)1便于 loss 函数加权。4.3 导出为 PyTorch 可加载的 HDF5兼顾速度与兼容性NetCDF 在 Python 中读取稍慢且部分嵌入式设备不支持。HDF5 是更通用的选择import h5py with h5py.File(china_temp_2020.h5, w) as f: f.create_dataset(lat, datagdf.geometry.y.values, dtypef4) f.create_dataset(lon, datagdf.geometry.x.values, dtypef4) f.create_dataset(temp_point, datagdf[temp_2020].values, dtypef4) f.create_dataset(temp_grid, datagdf[grid_temp].values, dtypef4) f.create_dataset(error, datagdf[abs_error].values, dtypef4) f.create_dataset(spatial_weight, datagdf[spatial_weight].values, dtypef4) # 添加 attributes 保留元信息 f[lat].attrs[units] degrees_north f[lon].attrs[units] degrees_east f[temp_point].attrs[units] degree_C # ... 其他变量同理PyTorch Dataset 示例class TempDataset(torch.utils.data.Dataset): def __init__(self, h5_path): self.f h5py.File(h5_path, r) self.keys [lat, lon, temp_point, temp_grid, error, spatial_weight] def __len__(self): return len(self.f[lat]) def __getitem__(self, idx): sample {} for k in self.keys: sample[k] torch.tensor(self.f[k][idx], dtypetorch.float32) return sample5. 避坑点面融合中 5 个血泪经验换来的高频翻车点这些不是教科书错误而是我在三个省级气候模型项目中亲手踩过的坑每一条都曾导致模型训练发散或论文被审稿人质疑数据可靠性。5.1 现象rasterio.sample()返回全nan原因栅格的nodata值未被正确识别。grid_0.1deg.tif的nodata可能是-9999、-32767或0但rasterio默认不读取sample时将 NoData 区域视为有效值而实际该位置无数据返回nan。解决务必在rasterio.open()后打印src.nodata并在reproject时显式传入dst_nodatasrc.nodata。若src.nodata is None用np.percentile(data, 0.1)估算下限值设为 nodata。5.2 现象点数据lon/lat与栅格bounds明显错位如点在 120°E栅格却显示 119.95°E原因grid_0.1deg.tif的transform存储的是左上角坐标但bounds是计算得出。某些国产栅格工具如 ArcGIS Export Raster会错误写入transform导致bounds与实际像元中心偏移半个像元。解决不用bounds用transform * (col0.5, row0.5)精确计算每个像元中心坐标。验证方法取栅格左上角像元计算其中心lon transform.c transform.a*0.5,lat transform.f transform.e*0.5与points.csv中最近站对比。5.3 现象Voronoi权重计算后青藏高原站点权重为 0原因Voronoi在凸包外生成无限区域intersection后面积为 0。shapely的polygon.area对退化多边形返回 0而非nan。解决在clipped.area后加判断if clipped.area 1e-6否则赋一个极小正值如1e-4避免除零。更鲁棒的做法是改用scikit-learn的NearestNeighbors计算 k5 邻居距离的倒数作为权重。5.4 现象NetCDF 中temp_point与temp_grid单位不一致一个 ℃一个 K原因部分再分析产品如 ERA5发布的是 Kelvin而points.csv是 ℃。ZIP 包未声明单位靠文件名“气温”想当然。解决永远用gdalinfo -stats grid_0.1deg.tif查看栅格统计值。若Min-273.15基本是 Kelvin若Min-50则是 ℃。转换公式K ℃ 273.15。在写入 NetCDF 前统一为 ℃并明确写入units属性。5.5 现象模型训练时 loss 突然爆炸debug 发现error变量含inf原因gdf[grid_temp]中有nanabs_error abs(a - b)在bnan时返回nan但nan参与 loss 计算会传播为inf。解决在计算abs_error前强制清洗gdf[grid_temp] gdf[grid_temp].fillna(gdf[grid_temp].median())或更合理——用gdf gdf.dropna(subset[grid_temp])剔除无法插值的站点并记录剔除数量通常 5% 属正常。6. 进阶技巧用空间残差图定位模型改进方向做完点面融合别急着扔进模型。真正的价值在于把error变量当诊断工具——它不是噪声而是揭示插值算法弱点的 X 光片。6.1 绘制全国空间残差分布图import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig plt.figure(figsize(12, 8)) ax plt.axes(projectionccrs.PlateCarree()) ax.set_extent([73, 135, 18, 54], crsccrs.PlateCarree()) # 绘制残差点图大小残差绝对值颜色正负 scatter ax.scatter( gdf.geometry.x, gdf.geometry.y, cgdf[abs_error], sgdf[abs_error] * 20, # 放大可视化效果 cmapRdBu_r, transformccrs.PlateCarree(), alpha0.7, edgecolorsblack, linewidth0.2 ) # 添加国界 ax.add_feature(cfeature.BORDERS, linestyle-, linewidth0.8) ax.add_feature(cfeature.COASTLINE, linewidth0.8) plt.colorbar(scatter, axax, label|Observed - Interpolated| (℃)) plt.title(2020 China Annual Mean Temperature Interpolation Error) plt.show()你会立刻发现误差高值2℃密集出现在三条带上——横断山脉、天山南麓、长白山东坡。这不是随机噪声而是地形强迫导致的插值失效区现有 0.1° 栅格无法解析 1km 级别的山谷风、逆温层。6.2 构建地形修正因子用 DEM 提升插值精度既然误差与地形强相关就把它变成特征# 下载 SRTM 30m DEM中国全境 # 使用 gdal.Warp 重采样到 0.1°与气温栅格同分辨率 # 此处假设已得 dem_0.1deg.tif with rasterio.open(dem_0.1deg.tif) as dem_src: # 提取每个站点的海拔 dem_samples list(sample_gen(dem_src, coords)) gdf[elevation] [s[0] for s in dem_samples] # 计算残差与海拔的线性关系 from sklearn.linear_model import LinearRegression X gdf[[elevation]].dropna() y gdf.loc[X.index, abs_error] model LinearRegression().fit(X, y) print(Elevation coefficient:, model.coef_[0]) # 通常为正表明海拔越高误差越大 # 生成地形修正栅格用 elevation 拟合残差再反推修正场 # 实际项目中用 GAM 或 RF 效果更好但 LinearRegression 已揭示核心规律6.3 误差驱动的模型架构选择何时该放弃 CNN如果你的下游任务是“预测未来某站气温”那么error分布告诉你在地形复杂区单纯的空间卷积CNN无法捕捉海拔-温度非线性关系。此时应切换为Graph Neural Network把气象站当作图节点用elevation、distance_to_coast、slope作为边权重——这正是 2023 年《Climate Dynamics》一篇高引论文的突破点。我最后养成的习惯是每次拿到新的气象数据 ZIP第一件事不是建模而是跑一遍points.csvgrid_xxx.tif的残差分析。那张红蓝斑驳的全国图比任何指标都诚实——它告诉我哪里该加物理约束哪里该换网络结构哪里该去野外补测。数据融合不是技术流程而是和地球对话的翻译过程。希望帮到你。本文还有配套的精品资源点击获取