简介中国2020年均气温数据点加栅格压缩包是一份面向ArcGIS用户的NCDC处理结果包含2020年全国年均气温站点观测点与栅格面两类数据。它既可以用于气象专题制图、区域气温空间插值也可以支撑城市热环境分析、农业气候区划等常见研究场景适合GIS学习者、地理信息相关专业学生以及需要全国尺度气温底图数据的科研人员使用。压缩包内共包含10个文件大小约100KB核心部分是tif格式的年均气温栅格和shp格式的站点点数据同时配有prj投影定义、tfw地理配准信息、sbn/sbx空间索引、ovr金字塔、dbf属性表等辅助文件解压后可在ArcGIS中直接加载坐标系显示正确浏览与编辑都非常高效。目前已有1636人学习下载具有一定的实用参考价值。站点点与栅格双格式并存既能实现点位数值查询又能通过栅格快速生成连续温度表面省去了自行下载、整理和预处理气象数据的步骤适合直接用于制图、插值对比或作为其他空间分析的基础数据。1. 中国2020年均气温数据点加栅格解压前先看清数据边界拿到一个名为「中国2020年均气温数据点加栅格.zip」的压缩包第一反应不应该是双击解压而是先确认里面装的到底是“点”还是“栅格”还是两者都有。气象业务里“数据点”通常是国家级或省级气象站点的观测值记录站点经纬度和2020年整年的平均气温而“栅格”是已经过空间插值或模型输出的连续温度面分辨率、投影、单位可能完全不一致。常见问题是直接把站点数据和已有栅格混在一起统计结果省界对不上、像元单位差十倍等发现时已经浪费了大半天。这个主题适合 GIS 数据处理、气候统计和环境建模的从业者核心思路是先把点状数据清洗成可插值的输入再生成栅格最后把栅格统一到同一坐标体系并裁剪归档。先花五分钟看文件清单比什么都重要。2. 拆 zip 与数据点清洗中文编码、时间字段、站点坐标2.1 先用 unzip -l 检查清单避免解压出一堆乱码文件zip 压缩包在跨平台传输时最容易出问题的不是数据本身而是文件名的编码。Windows 资源管理器生成的 zip 常用 GBK 编码文件名Linux 和 macOS 默认按 UTF-8 解码所以直接在 Linux 下unzip经常会看到一堆中文乱码文件名。更麻烦的是有些包在压缩时把目录结构压坏了直接解压会得到几十个文件平铺在当前目录。我一般会先做一次“只读不解压”的检查unzip -l 中国2020年均气温数据点加栅格.zip | head -30-l是 list 的缩写只列出压缩包中的文件名、原始大小、压缩后大小和修改时间不实际写入磁盘。head -30用来控制输出行数避免包内文件过多时终端刷屏。看到清单后再判断是否需要全部解压还是只抽出需要的 CSV 和 TIF。如果文件名是乱码先试unzip -O gbk -l这是 Info-ZIP 对非 UTF-8 编码的兼容选项。部分发行版编译时没有开启该选项那就不必纠结直接用 7-Zip 解压或者写个 Python 脚本按字节编码来重命名。若提示error read zip archive说明压缩包损坏或使用了非常规加密头先跑一遍zip -T 中国2020年均气温数据点加栅格.zip验证完整性任何所谓“密码破解工具”都不该成为首选先向数据来源方要完整文件。2.2 数据点表格先看最小值、最大值和单位气温数据最常见的坑不是读不出来而是单位没换算。气象观测的原始数据中0.1℃甚至 0.01℃ 都是常见存储单位列名可能叫TEMP、TAVG、MEAN_TEMP不打开看绝对值根本不知道真实含义。先用 pandas 读取并做快速描述性统计import pandas as pd df pd.read_csv(points_2020.csv, encodingutf-8) print(df.head()) print(df[[longitude, latitude, temp]].describe()) if df[temp].abs().max() 80: df[temp] df[temp] / 10 print(temp units adjusted from 0.1°C to °C)代码先按 UTF-8 读文件如果报UnicodeDecodeError就改成encodinggbk或encodingutf-8-sig因为部分 Windows 导出的 CSV 带 BOM 头。describe()输出均值、标准差、最小值和最大值看到极端值之后再决定要不要做单位换算。if df[temp].abs().max() 80是一个快速判据中国 2020 年年平均气温不会超过 30℃ 太多如果绝对值大于 80基本可以断定是 0.1℃ 精度。写完单位换算后接着按站号或经纬度去重并去掉坐标或温度为空的行df df.drop_duplicates(subset[station_id]) df df.dropna(subset[longitude, latitude, temp])station_id是站点标识列如果源数据没有这列就改成subset[longitude, latitude]防止同一站点的重复观测被当成多个样本参与插值。2.3 坐标参考判断数据点是经纬度还是投影坐标有些 csv 里的列名是x和y不代表它就是投影坐标有些则明明是经纬度却叫x/y。判断方法很简单看数值范围。中国区域的经纬度大致落在经度 73°E135°E、纬度 18°N53°N如果是投影坐标X 可能是百公里到千公里级别的米制值。这里列一张快速判断表字段取值特征可能坐标系后续处理建议longitude 在 73135latitude 在 1853WGS84 或 CGCS2000 地理坐标系先按 EPSG:4326 或 EPSG:4490 处理x 在 500000几百万y 在几百万上千万高斯-克吕格投影需要按中央经线或带号找到对应 EPSGx 在几十万以内y 在几十万以内可能是 Web Mercator 或地方坐标系谨慎先查看 .prj 文件如果数据点是 Shapefile可以直接用ogrinfo看它的坐标系声明ogrinfo -al -so points.shp-so是 summary only只输出图层摘要不会逐条打印要素属性。输出中会有一行Layer SRS WKT从这里能看到完整的坐标参考描述。如果没有 .prj 文件就按上表结合数值范围推测。这个步骤不能省因为后面所有插值和栅格转换都建立在“我知道当前坐标是什么”的基础上。3. 从数据点插值生成平均气温栅格IDW 到克里金的参数取舍3.1 网格分辨率先从 0.1° 起步站点气温插值的第一步不是选算法而是先定输出栅格的像元大小。栅格分辨率越低运算越快但会掩盖局部地形对气温的影响分辨率越高插值器越容易在站点稀疏区制造虚假细节。我一般先用 0.1°×0.1° 做一版全国试验大约相当于 10km 格距检查误差和空间连续性后再考虑加密到 0.05°。建立经纬度网格import numpy as np lon_min, lon_max, lat_min, lat_max 73, 135, 18, 54 res 0.1 lon_grid np.arange(lon_min, lon_max, res) lat_grid np.arange(lat_min, lat_max, res)这里用了最简单的地理坐标系网格单位是度。np.arange的终点是不包含边界所以实际网格到不了 135°E 和 54°N后续可以用栅格的范围重新补齐。因为中国国界不是规则的矩形这里先建一个规则网格最后再裁剪。如果要更接近实际面积可以换用 1km 的投影坐标系网格但在全国尺度下0.1° 已经能反映大尺度温度格局。3.2 IDW 与克里金怎么选变差函数、搜索半径、各向异性插值方法的选择取决于站点密度和地形复杂度。IDW 形式简单适合快速查看空间走势克里金会给出一个伴随的误差面更适合要求结果中带诊断信息的场景。ArcGIS Pro 里常用的 IDW、克里金以及“地形转栅格”在这个问题上不是互相替代的关系而是根据场景取舍。方法适用场景核心参数常见坑IDW站点密、地形平坦power2搜索半径 1°2°功率过大时站点周围出现“牛眼”普通克里金站点较均匀需要误差估计variogram_modelsphericalnlags6站点少于 50 个时拟合极不稳定样条函数追求平滑连续表面tension 控制弯曲程度插值结果可能超过观测值范围地形转栅格有高程栅格辅助建模drainage enforcement 需关闭针对地形水文设计气温插值要慎用Python 里可以用pykrige做普通克里金from pykrige.ok import OrdinaryKriging ok OrdinaryKriging( lonsdf[longitude].values, latsdf[latitude].values, datadf[temp].values, variogram_modelspherical, nlags6, weightTrue, ) z, ss ok.execute(grid, lon_grid, lat_grid)参数里variogram_modelspherical是最常用的球状变差函数模型适合气温这种空间自相关性随着距离衰减平稳的数据nlags6表示把站点两两之间的距离分成 6 组来统计经验变差函数组数太少会掩盖近距离的空间结构组数太多会让每个组的样本不足。weightTrue让变差函数拟合时按点对数加权。execute(grid, ...)返回两个数组z是插值结果ss是克里金方差。注意z的形状是(len(lat_grid), len(lon_grid))写 GeoTIFF 或绘图时不要转错转置。3.3 站点密度差异带来的空值问题中国站点分布极不均匀东部几十公里就有一个站青藏高原和新疆腹地可能几百公里都没有一个。插值器不会主动告诉你某个位置离最近站点已经超过合理范围它只会把远处的外推值照样算出来甚至给出一个很平滑但完全不可信的温度。所以插值后必须做“近距离约束”from scipy.spatial import cKDTree coords df[[longitude, latitude]].values grid_coords np.column_stack([lon_grid.ravel(), lat_grid.ravel()]) tree cKDTree(coords) dist, _ tree.query(grid_coords, k1) dist dist.reshape(z.shape) z np.where(dist 1.0, z, np.nan)思路是用 KDTree 找到每个网格点的最近站点距离然后把距离超过阈值的插值结果置为 NaN。阈值 1.0 表示 1 度左右如果只关心东部地区可以缩小到 0.5如果是全国尺度且要保留西部稀疏区的大致形态可以放大到 2.0。这个距离筛选比任何插值算法调参都重要它把“不可知”的区域诚实地标记为空值而不是让读者误以为那里有可靠数据。4. 栅格后处理WGS84 转 CGCS2000、裁剪与 ArcGIS Pro 栅格合并4.1 统一投影和像元gdalwarp 的常用参数不同来源的栅格产品坐标系可能是 WGS84也可能是 CGCS2000甚至同一批分块栅格里混着两套坐标。把数据点插值得到的栅格转成 CGCS2000 时最直接的做法是gdalwarp。这里涉及到一个常见长尾问题栅格影像 WGS84 坐标系转 CGCS2000。gdalwarp -s_srs EPSG:4326 -t_srs EPSG:4490 \ -tr 0.1 0.1 -r bilinear -overwrite \ temp_2020_wgs84.tif temp_2020_cgcs2000.tif-s_srs EPSG:4326指定输入数据是 WGS84 地理坐标系-t_srs EPSG:4490指定输出为 CGCS2000 地理坐标系。-tr 0.1 0.1是目标像元尺寸单位与目标坐标系一致投影后仍是经纬度。-r bilinear选择双线性重采样适合连续型栅格如果做土地利用分类就不该用它应该用最近邻。-overwrite表示允许覆盖已存在的输出文件。有个细节容易忽略WGS84 与 CGCS2000 在大部分地区差异只有米级但如果源数据是投影坐标直接按经纬度转换会差得离谱。所以先确认源数据的坐标参考再决定是否需要-s_srs。4.2 用省界或国界裁剪栅格拿到全国温度栅格后第一版范围是一个矩形落到边界外全是空值或插值外推值。这时候需要用中国边界矢量做裁剪gdalwarp -s_srs EPSG:4490 -t_srs EPSG:4490 \ -cutline 全国边界.shp -crop_to_cutline \ temp_2020_cgcs2000.tif temp_2020_clip.tif-cutline后面接边界矢量文件-crop_to_cutline让输出范围严格等于边界范围而不是边界所在矩形的四至。注意裁剪前要确认边界矢量和栅格坐标一致否则 gdalwarp 会自动做一次坐标系变换边界越复杂重投影耗时就越高。裁剪完用gdalinfo temp_2020_clip.tif查看Size和Corner Coordinates检查范围是否合理。4.3 ArcGIS Pro 栅格合并与 NoData 处理如果产品本身是分省或分块栅格后面就要做合并。ArcGIS Pro 里叫“Mosaic to New Raster”也就是常说的 ArcGIS Pro 栅格合并GDAL 里对应的是gdalbuildvrt加gdal_translate。gdalbuildvrt temp_2020.vrt 分省/*.tif gdal_translate -co COMPRESSDEFLATE temp_2020.vrt temp_2020_all.tifgdalbuildvrt不复制实际像素只是把多个分块文件的引用写进一个 VRT 虚拟目录所以速度很快。gdal_translate再把虚拟目录落成一个正式 GeoTIFF-co COMPRESSDEFLATE指定 LZ77 压缩能显著减小 zip 文件体积。合并前必须检查各分块的 NoData 值是否一致如果一个文件用-9999另一个用0合并后就会出现大片黑色或白色伪值。操作工具关键设置注意事项分块合并gdalbuildvrt不设-r时默认取第一个块各块投影必须一致镶嵌导出gdal_translate-co COMPRESSDEFLATENoData 统一图形界面合并ArcGIS Pro Mosaic to New RasterPixel Type、NoData Value、Mosaic Operator重叠区选 BLEND 或 FIRST合并后检查gdalinfo-stats计算统计值观察 min/max 是否异常gdalbuildvrt拼接时如果块与块之间有重叠默认的行为可能会生硬这时需要回到 ArcGIS Pro 中设置 Mosaic Operator。对气温这种渐变栅格BLEND比FIRST更能消除接缝。5. 精度验证与可视化残差、RMSE 与专题图5.1 留一交叉验证把每根数据点当作隐藏测试插值结果不能只靠“看起来像不像”需要用站点数据本身做留一交叉验证。所谓留一法就是从全部站点中依次剔除一个点用其余点点位重新插值再预测被剔除点的温度最后统计预测与实测的误差。这个过程能直接回答“如果某个站点没有观测插值能猜得多准”。import numpy as np from pykrige.ok import OrdinaryKriging lon df[longitude].values lat df[latitude].values t df[temp].values err [] for i in range(len(t)): mask np.ones(len(t), dtypebool) mask[i] False ok OrdinaryKriging(lon[mask], lat[mask], t[mask], variogram_modelspherical, nlags5) pred, _ ok.execute(points, lon[i], lat[i]) err.append(t[i] - float(pred[0])) err np.array(err) rmse np.sqrt(np.mean(err**2)) mae np.mean(np.abs(err)) print(RMSE , rmse, MAE , mae)这个循环的代价很高几千个站点就相当于几千次克里金拟合运行时间可能从几分钟到数小时。所以实际项目中如果站点超过 300 个我会改成 5 折分组验证也就是随机分成五组每次拿其中一组作为验证集剩下四组用于插值。err数组还能进一步分析空间分布比如把残差按经纬度画成散点图能看出西部和东部哪个区域误差大。RMSE 对异常值敏感MAE 更稳健两者一起报告就够了。5.2 用 rasterio 读取栅格并绘制平均气温专题图验证完再出图。用rasterio读裁剪后的栅格matplotlib出图是最常见的一种做法import rasterio from matplotlib import pyplot as plt src rasterio.open(temp_2020_clip.tif) temp src.read(1, maskedTrue) fig, ax plt.subplots(figsize(10, 8)) im ax.imshow(temp, cmapSpectral_r, originupper) contour ax.contour(temp, levels[0, 10, 15, 20], colorsblack, linewidths0.5) plt.colorbar(im, label2020 mean temperature (°C)) plt.savefig(temp_2020_map.png, dpi200)maskedTrue会把 NoData 自动转成 NumPy 的掩膜数组绘图时不会呈现为 0 值。originupper与 GeoTIFF 的行列方向匹配因为很多栅格数据的第一行是北侧。等值线叠加时可以限定在特定温度范围比如把 0℃、10℃、15℃ 和 20℃ 线标出来方便判断积温和熟制带。如果要在图上加国界可以引入 Cartopy 的 shapefile 支持但单独输出一个纯净的栅格图通常更安全。5.3 和已有的官方栅格产品对比偏差如果压缩包里原本就带一个官方平均气温栅格不要急着扔掉用它对插值结果做差值检验。GDAL 的栅格计算器可以一行命令完成gdal_calc.py -A temp_2020_clip.tif -B official_2020.tif \ --outfilediff.tif --calcA-B --NoDataValue-9999gdal_calc.py是 GDAL 自带的 Python 脚本-A和-B分别给两个输入栅格命名--calcA-B支持简单的 NumPy 表达式--NoDataValue-9999保证差值图在空值区域也不会污染统计。差值图如果出现大范围正偏差说明插值结果整体偏高先回头检查站点温度单位如果正负偏差交错但幅度大于 2℃再考虑是插值算法问题还是数据年份口径不同。差值图的方差比 RMSE 更能表达空间结构。6. 重新封 zip目录约定、坐标系文件和元数据6.1 目录约定避免第三个人接手时再猜一遍处理完的数据如果还是以一个 zip 压缩包发出去那么包内目录一定要和原始包有区别。我一般会按“数据点、栅格、文档、元数据”四类组织import zipfile with zipfile.ZipFile(中国2020年均气温数据点加栅格_processed.zip, w, compressionzipfile.ZIP_DEFLATED, compresslevel9) as zf: zf.write(points_2020.csv, data/points/points_2020.csv) zf.write(temp_2020_clip.tif, data/grids/temp_2020_clip.tif) zf.write(README.txt, doc/README.txt) zf.write(metadata.xml, metadata.xml)ZIP_DEFLATED表示普通无损压缩compresslevel9是最高压缩级别。GeoTIFF 本身已经压缩过再压不会小多少但对于 CSV 和文档来说收益明显。更关键的是文件名路径CSV、TIF、README 分别放在不同目录避免别人在解压时把所有文件挤在一个目录下面。归档时不要漏掉坐标系说明。如果数据点是 Shapefile.prj文件必须和.shp放在一起如果只是 CSV就在 README.txt 里写清楚“经度纬度坐标基于 EPSG:4490温度单位 ℃”。6.2 验证归档不是“能打开”就叫完整压缩包写完最后做一次完整性校验7z t 中国2020年均气温数据点加栅格_processed.zip7z t会逐个文件解压到内存并比对 CRC 校验值输出 “Everything is Ok” 才表示文件没有损坏。但这也只验证了 zip 结构完整不代表栅格和坐标系没问题。更严格的做法是用 GDAL 的 /vsizip/ 虚拟文件系统直接读取包内 TIFgdalinfo /vsizip/中国2020年均气温数据点加栅格_processed.zip/data/grids/temp_2020_clip.tif/vsizip/是 GDAL 内置的虚拟文件系统能不解压直接读取 zip 内的栅格或矢量文件。如果这行命令能正确打印出投影、像元大小和 NoData 值说明归档后的成果在 GIS 软件里也能顺利打开。顺手把zip -sf列出压缩包内文件的输出存进验证记录整个数据点加栅格的处理流程才算收尾。本文还有配套的精品资源点击获取