简介成渝城市群1980—2015年逐年气温栅格数据面向自然地理、人文地理、环境科学、生态学等专业的学生和科研人员可为区域气候演变、城市群热环境效应及生态过程模拟提供基础数据支撑。数据集源自全国2400多个气象站点的逐日观测记录经整理、计算与空间插值后生成并完成研究区裁剪与数值换算处理栅格值扩大10倍除以10即为实际摄氏度。数据时间跨度为1980年至2015年涵盖逐年平均气温与降水量指标用户可在ArcGIS、QGIS等软件中直接加载调用。资源共114个文件包含36个tif栅格及配套tfw坐标文件、xml元数据另有少量dat、nit等辅助文件压缩包大小约21.49MB目录结构简洁便于按年份检索使用。目前已有622人学习下载尤其适合需要长时序区域气象栅格数据开展GIS分析与遥感应用的研究者可直接用于趋势分析、制图或模型输入为成渝地区气候研究提供统一投影坐标的栅格产品省去站点数据整理与插值步骤。1. 拿到成渝城市群气温栅格先别急着双击 .dat第一次拿到这份数据的人大概率会被arc0000.dat、arc.dir、arc0000.nit、Tem1991.tfw这一串文件搞蒙。这不是简单的二进制文件而是 ArcGIS 传统的 Grid 栅格格式一个文件夹就是一个完整的栅格数据集。数据覆盖成渝城市群时间跨度为 1980 到 2015 年包含逐年平均气温和降水量两种要素属性值已经放大 10 倍除以 10 才是真实的摄氏度或毫米数。对自然地理、生态和环境专业的研究者来说这套数据可以直接用于城市热岛、气候带迁移、农业热量资源分析等场景前提是你要知道怎么把它正确读进来、换算、裁剪再做时空趋势分析。下面从文件结构开始一步一步拆解。2. 从 arc0000.dat 到可计算的温度场栅格格式解析与 GDAL 读取2.1 ArcGIS Grid 的文件构成与识别要点ArcGIS Grid 是一种老牌栅格格式分为二进制 Grid 和 ASCII Grid 两种。这份数据里的arc0000.dat是二进制 Grid 的核心数据文件存的是像元值矩阵arc.dir是索引文件记录栅格的行列数、像元大小、空间范围arc0000.nit是头信息文件包含 NoData 值、数据精度等元数据log是历史日志。带有.tfw的文件则是世界文件例如Tem1991.tfw用于把像元行列坐标转换成真实地理坐标实际上这部分信息在 Grid 的arc.dir中已经存在.tfw是给其他软件做外部配准用的。判断一个目录是不是 ArcGIS Grid最重要的标志就是看它里面是否有arc0000.dat和arc.dir。在 Windows 上ArcGIS 会把这个目录显示为一个栅格图层但如果你直接打开该文件夹里的单个.dat文件用 Notepad 看会得到乱码。正确做法是用 GDAL 或 ArcPy 打开整个目录路径而不是打开具体文件。在读取之前先确认数据的空间范围与投影。如果你有 ArcGIS 或者 QGIS直接在 QGIS 中添加该目录QGIS 会通过 GDAL 驱动识别它。但为了后续批量处理我还是建议用 Python 方式操作这也是多数气候数据处理流程中最省事的一环。2.2 用 Python GDAL 把温度读成 ndarrayGDAL 的 Grid 驱动可以自动识别这种目录式数据。假设你拥有一个Tem1991目录里面放着上述文件用下面的代码就能读到矩阵from osgeo import gdal import numpy as np path Tem1991 ds gdal.Open(path, gdal.GA_ReadOnly) if ds is None: raise IOError(无法打开栅格数据检查目录结构是否完整) print(行列数:, ds.RasterXSize, ds.RasterYSize) print(波段数:, ds.RasterCount) print(投影:, ds.GetProjection()) print(仿射变换参数:, ds.GetGeoTransform()) band ds.GetRasterBand(1) arr band.ReadAsArray() print(原始数组dtype:, arr.dtype, 最小值:, arr.min(), 最大值:, arr.max())这段代码的核心是gdal.Open注意传入的路径是Tem1991这个目录而不是arc0000.dat。如果路径正确GDAL 会自动识别 Grid 驱动。GetGeoTransform()返回六个值的元组分别是左上角 X 坐标、像元宽度、旋转项通常为 0、左上角 Y 坐标、旋转项、像元高度通常为负值。有了这组参数之后就可以把像元行列号转换成经纬度或投影坐标。很多人在这一步会直接读取.dat文件导致失败或者得到一维二进制流所以一定要从目录层面打开。2.3 数值还原与 NoData 处理数据说明里明确指出所有值都扩大了 10 倍除以 10 才是实际温度。因此读取后必须做一次换算# 将原始值转为浮点型并统一除以10得到摄氏度 temp_celsius arr.astype(np.float32) / 10.0 print(温度范围℃:, np.nanmin(temp_celsius), np.nanmax(temp_celsius))这里有一个容易被忽略的问题原始数组如果是整型除以 10 后得到的是小数所以要先astype(np.float32)再做除法否则整数除以整数会丢失小数位。另一个关键是 NoData 值。有些栅格在无数据区域使用-9999或0表示缺测如果不处理直接除以 10 后会得到无效的负值或 0干扰后续统计。建议在读取波段后立即检查 NoData 值nodata band.GetNoDataValue() print(NoData 原始值:, nodata) if nodata is not None: temp_celsius[temp_celsius nodata / 10.0] np.nan注意这里用nodata / 10.0与换算后的数组比较因为原始数组还没换算时NoData 值是原始编码。如果在换算前标记就直接比较原始数组然后赋 NaN 再换算。养成这个习惯后续做区域均值和趋势分析时才不会出现异常低值。降水数据同理只是单位从摄氏度变成毫米同样也要注意放大 10 倍的问题。3. 裁剪、重投影与像元归一把原始栅格变成分析可用数据3.1 为什么需要裁剪和重投影成渝城市群是一个行政区划概念而这份栅格数据很可能覆盖的是整个西南地区甚至全国。直接使用全区域计算不仅效率低而且边界外像元会干扰城市群尺度的统计。裁剪势在必行。另一个问题是投影。从数据产出来看原始数据很可能采用了 Albers 等积投影或兰伯特投影以便保证面积不变。但如果你之后要跟站点数据、气象插值数据合并或者用seaborn、geopandas画图通常需要把投影统一到 WGS84 地理坐标系EPSG:4326上这样经纬度才是直观的坐标。需要注意重投影会改变像元值吗如果使用最近邻采样像元值不会插值改变但面积和像元大小会变化如果使用双线性插值则会对温度值做平滑这取决于你的用途。气候栅格分析一般都建议用最近邻或者先完成数值统计再做投影转换。3.2 用 Rasterio 按成渝城市群边界裁剪Rasterio 是 GDAL 的 Python 封装处理栅格裁剪比直接写 GDAL 更简洁。先备好一份成渝城市群的矢量边界比如chengyu_boundary.shp然后用rasterio.mask.mask完成裁剪import geopandas as gpd import rasterio from rasterio.mask import mask as rio_mask boundary gpd.read_file(chengyu_boundary.shp) # 确保边界文件与栅格的CRS一致不一致时先转换 boundary boundary.to_crs(EPSG:4326) with rasterio.open(Tem1991) as src: # 读取原始投影若边界与栅格投影不同此处需要先把边界转到栅格投影 geoms boundary.geometry.values out_image, out_transform rio_mask(src, geoms, cropTrue, nodatasrc.nodata) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(Tem1991_chengyu.tif, w, **out_meta) as dst: dst.write(out_image)代码逻辑分三步第一步读取矢量边界第二步打开栅格并调用mask第三步把裁剪结果写为 GeoTIFF。cropTrue表示让输出范围严格贴合边界外接矩形nodata会继承原始栅格的 NoData 值。要注意的是mask函数要求边界的坐标系与栅格一致所以最好先把边界转换到栅格自身坐标系或者用rasterio.warp.transform_geom做几何转换。很多人在这一步报错就是因为投影不一致。另外如果你希望裁剪后的栅格保留原始投影而不是边界投影就不要提前把边界转成 EPSG:4326而是让 GDAL 自己处理。3.3 栅格统计区域均值怎么算才对裁剪完成后计算成渝城市群范围内的年均温均值最简单的办法是把像元值取平均但不要直接对裁剪后的数组做np.mean因为边界外的 NoData 会变成 0 或者无效值。正确姿势是使用np.nanmean并结合 NoData 掩膜。import numpy as np with rasterio.open(Tem1991_chengyu.tif) as src: band src.read(1, maskedTrue) # maskedTrue 会把NoData自动转为masked # band现在是一个masked array计算平均值时自动忽略无效值 region_mean band.mean() print(1991年成渝城市群平均气温由放大10倍数值还原:, float(region_mean) / 10.0, ℃)read(1, maskedTrue)返回一个numpy.ma.MaskedArray所有 NoData 都被隐藏在 mask 中mean()会跳过它们。如果不用 masked 模式就需要自己构建掩膜data src.read(1) nodata src.nodata valid data ! nodata mean_val data[valid].astype(np.float32).mean() / 10.0这两种方式都可以但maskedTrue在大规模计算时更直观。计算区域均值时还有一个容易被忽视的因素像元面积。如果栅格是等经纬度投影高纬度像元对应的面积比低纬度小此时计算区域平均温度应当做面积加权。但在成渝城市群这样一个南北跨度约 3 到 4 度的范围内面积差异很小做不等权校正的影响极其有限。如果你处理的是全国尺度数据建议用xarray的cos(latitude)加权。这里先按普通均值处理后续趋势分析也沿用这一口径确保不同年份之间统计方法一致。4. 1980—2015 年气温时空趋势斜率、显著性检验与可视化4.1 像元级线性趋势计算拿到 36 年年均温栅格后最想回答的问题通常是成渝城市群哪里在升温升温最剧烈的地方在哪这就需要对每个像元做时间序列回归。假设你已经把 1980 到 2015 年全部裁剪后的栅格读入一个三维数组形态为(36, rows, cols)那么逐像元的线性趋势可以用最小二乘拟合来求。为了效率不要用 Python 双重循环遍历所有像元而是把空间维展开成二维一次性计算years np.arange(1980, 2016, dtypenp.float64) # stack 形状:(36, rows, cols) # 转成 (36, rows*cols) n_years, rows, cols stack.shape flat stack.reshape(n_years, -1) # (36, N) # 中心化年份减少数值计算误差 years_c years - years.mean() denom (years_c ** 2).sum() # 计算每个像元的斜率 # 斜率 Σ( (year - year_mean) * (value - value_mean) ) / Σ(year - year_mean)^2 # 先计算每个像元在时间维上的方差 slope np.full(rows*cols, np.nan, dtypenp.float64) valid_mask np.all(~np.isnan(flat), axis0) if valid_mask.sum() 0: xn years_c.reshape(-1, 1) - years_c.mean()上面的写法有点绕更直接的方式是用向量化矩阵乘法# 标准化年份向量 x years - years.mean() x_std x / np.sqrt((x ** 2).sum()) # 对时间维中心化后的数据做回归 flat_centered flat - np.nanmean(flat, axis0, keepdimsTrue) # 斜率 sum(x_std * centered) / sqrt(sum(x^2)) , 由于x_std已经归一化 # 实际手动实现 y_demeaned flat - flat.mean(axis0, keepdimsTrue) # (36,N) slope_vec (x_std.reshape(-1,1) * y_demeaned).sum(axis0) / np.sqrt((x ** 2).sum())最终slope_vec就是每个像元的气温年倾率单位为℃/年。注意这里的数据必须已经是真温度也就是除以 10 后的数值。如果直接用原始编码值得到的斜率是 0.2 ℃/年 的 10 倍容易误判。计算完成后把slope_vecreshape 回二维就得到升温速率空间分布图。这个斜率就是经典的 Sen 斜率的一种简单线性估计适用没有异常值的年度序列。如果数据中有个别年份缺测np.isnan的掩膜会过滤掉但需要注意至少要有 30 年以上有效数据否则拟合结果稳健性差。4.2 区域平均温度序列的趋势检验像元级趋势能告诉你空间分布但政策或研究报告通常需要一条成渝城市群年平均气温区域序列并判断其变化是否统计显著。这里推荐 Mann-Kendall 趋势检验它不要求数据正态分布对异常值稳健。Python 中可用pymannkendall库实现pip install pymannkendall然后计算区域平均序列import pymannkendall as mk # 对每年stack取有效像元均值 annual_mean [] for i in range(stack.shape[0]): valid stack[i] annual_mean.append(np.nanmean(valid)) annual_mean np.array(annual_mean) # 单位℃ result mk.original_test(annual_mean, alpha0.05) print(MK趋势统计量:, result.Tau) print(p值:, result.p) print(趋势方向:, result.trend) print(斜率估计:, result.slope)pymannkendall的original_test返回一个结果对象其中trend会给出 increasing、decreasing 或 no trendslope则是 Sen 斜率单位是℃/年。相比普通最小二乘回归MK 检验不依赖线性假设能更保守地判断是否真的存在单调趋势。实测中成渝城市群年均温序列在 1980—2015 年间多呈现弱升温趋势但不同城市群子区域差异很大。如果你想分析四季或月平均气温趋势需要先把原始数据按季度合成再重复上述流程。4.3 趋势图与结果导出趋势结果最好用图来呈现。用matplotlib画两个图左图是像元趋势空间分布右图是区域平均时间序列与拟合线。先用basemap或cartopy绘制底图没有现成的城市群边界时直接使用裁剪后栅格的经纬度经纬度范围。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature lon np.linspace(transform[0], transform[0] transform[1]*cols, cols) lat np.linspace(transform[3] transform[5]*rows, transform[3], rows) fig plt.figure(figsize(12,5)) ax1 fig.add_subplot(1,2,1, projectionccrs.PlateCarree()) im ax1.imshow(slope_2d, extent[lon[0], lon[-1], lat[-1], lat[0]], cmapRdBu_r, transformccrs.PlateCarree()) ax1.add_feature(cfeature.COASTLINE) ax1.coastlines() plt.colorbar(im, axax1, shrink0.6, label℃/年) ax2 fig.add_subplot(1,2,2) ax2.scatter(years, annual_mean, s20, ck, alpha0.6) ax2.plot(years, slope* years intercept, r-, label线性拟合) ax2.set_xlabel(年份); ax2.set_ylabel(年均温℃) plt.show()绘制趋势图时要注意cartopy的imshow坐标参数extent需要按照左下角、右下角、右上角、左上角的经纬度顺序传入否则图像会上下颠倒。lat[-1]到lat[0]的写法是为了适应纬度递减的情况。趋势斜率图建议使用RdBu_r调色板红色代表升温蓝色代表降温符合气候报告的习惯。导出结果时可以直接用numpy.savetxt保存趋势矩阵或者用rasterio写成 GeoTIFF这样后续在 ArcGIS 中可以直接叠加查看。5. 批量转 NetCDF把 36 年温度数据装进一个文件5.1 用 xarray 合并逐年栅格逐年处理 36 个 GeoTIFF 非常麻烦而且每次读取都要重新打开文件。最理想的方式是把所有年份打包成一个多维 NetCDF维度为(time, latitude, longitude)之后用xarray一个文件搞定。首先要准备一个文件列表按年份排序然后逐个读取并构建坐标。import xarray as xr import numpy as np from osgeo import gdal file_list [fTem{year}_chengyu.tif for year in range(1980, 2016)] stack_list [] for year, path in zip(range(1980, 2016), file_list): ds gdal.Open(path) arr ds.ReadAsArray().astype(np.float32) / 10.0 gt ds.GetGeoTransform() # 计算每个像元的中心坐标注意gt[5]通常为负 lon gt[0] np.arange(ds.RasterXSize) * gt[1] gt[1] / 2 lat gt[3] np.arange(ds.RasterYSize) * gt[5] gt[5] / 2 stack_list.append(arr) ds None data np.array(stack_list) # (36, rows, cols)这里要注意gt[5]是负值表示纬度从北向南递减所以计算出的lat数组是递减的。xarray要求坐标单调递增因此需要翻转纬度轴data data[:, ::-1, :] # 翻转lat维度, 使其从小到大 lat lat[::-1] da xr.DataArray( data, dims(time, latitude, longitude), coords{time: list(range(1980, 2016)), latitude: lat, longitude: lon}, nametemperature, ) da.attrs[units] degC da.attrs[scale] original_values_divided_by_10xarray.DataArray在构建时如果坐标是浮点且不严格单调会自动报错。翻转后latitude变成从小到大的递增序列符合 CF 约定。如果你同时处理降水可以增加一个precipitation变量或者使用xr.Dataset包含两个变量。5.2 一次写出的 NetCDF 怎么用da.to_netcdf(chengyu_t2m_1980_2015.nc, enginenetcdf4)这行代码就会生成一个约几十兆的 NetCDF 文件包含全部 36 年的温度数据。之后再做任何分析都不需要再碰原始目录ds xr.open_dataset(chengyu_t2m_1980_2015.nc) temp ds[temperature] # (time, lat, lon) # 快速计算2000年之后的平均温度 recent_mean temp.sel(timeslice(2000, 2015)).mean(dimtime) # 快速提取某个点的温度序列 site_series temp.sel(longitude104.06, latitude30.67, methodnearest)xarray的切片、聚合、重采样能力比裸 NumPy 高效得多尤其是做多年平均、季节平均时直接用.groupby(time.year)或.resample(timeYS)即可。另外NetCDF 文件自带自描述信息在 Python 中读取不需要关心 NoDataxarray会将缺测值表示为NaN统计时会自动跳过。如果你要把数据交给 R 的ncdf4包处理或 Fortran 程序读取NetCDF 也是通用性最好的载体。这一步骤完成后原始数据的 10 倍放大问题已经彻底抹平后续所有分析的输入都是标准的物理量。本文还有配套的精品资源点击获取 成渝城市群气温栅格数据处理全流程：GDAL读取、裁剪与趋势分析
简介成渝城市群1980—2015年逐年气温栅格数据面向自然地理、人文地理、环境科学、生态学等专业的学生和科研人员可为区域气候演变、城市群热环境效应及生态过程模拟提供基础数据支撑。数据集源自全国2400多个气象站点的逐日观测记录经整理、计算与空间插值后生成并完成研究区裁剪与数值换算处理栅格值扩大10倍除以10即为实际摄氏度。数据时间跨度为1980年至2015年涵盖逐年平均气温与降水量指标用户可在ArcGIS、QGIS等软件中直接加载调用。资源共114个文件包含36个tif栅格及配套tfw坐标文件、xml元数据另有少量dat、nit等辅助文件压缩包大小约21.49MB目录结构简洁便于按年份检索使用。目前已有622人学习下载尤其适合需要长时序区域气象栅格数据开展GIS分析与遥感应用的研究者可直接用于趋势分析、制图或模型输入为成渝地区气候研究提供统一投影坐标的栅格产品省去站点数据整理与插值步骤。1. 拿到成渝城市群气温栅格先别急着双击 .dat第一次拿到这份数据的人大概率会被arc0000.dat、arc.dir、arc0000.nit、Tem1991.tfw这一串文件搞蒙。这不是简单的二进制文件而是 ArcGIS 传统的 Grid 栅格格式一个文件夹就是一个完整的栅格数据集。数据覆盖成渝城市群时间跨度为 1980 到 2015 年包含逐年平均气温和降水量两种要素属性值已经放大 10 倍除以 10 才是真实的摄氏度或毫米数。对自然地理、生态和环境专业的研究者来说这套数据可以直接用于城市热岛、气候带迁移、农业热量资源分析等场景前提是你要知道怎么把它正确读进来、换算、裁剪再做时空趋势分析。下面从文件结构开始一步一步拆解。2. 从 arc0000.dat 到可计算的温度场栅格格式解析与 GDAL 读取2.1 ArcGIS Grid 的文件构成与识别要点ArcGIS Grid 是一种老牌栅格格式分为二进制 Grid 和 ASCII Grid 两种。这份数据里的arc0000.dat是二进制 Grid 的核心数据文件存的是像元值矩阵arc.dir是索引文件记录栅格的行列数、像元大小、空间范围arc0000.nit是头信息文件包含 NoData 值、数据精度等元数据log是历史日志。带有.tfw的文件则是世界文件例如Tem1991.tfw用于把像元行列坐标转换成真实地理坐标实际上这部分信息在 Grid 的arc.dir中已经存在.tfw是给其他软件做外部配准用的。判断一个目录是不是 ArcGIS Grid最重要的标志就是看它里面是否有arc0000.dat和arc.dir。在 Windows 上ArcGIS 会把这个目录显示为一个栅格图层但如果你直接打开该文件夹里的单个.dat文件用 Notepad 看会得到乱码。正确做法是用 GDAL 或 ArcPy 打开整个目录路径而不是打开具体文件。在读取之前先确认数据的空间范围与投影。如果你有 ArcGIS 或者 QGIS直接在 QGIS 中添加该目录QGIS 会通过 GDAL 驱动识别它。但为了后续批量处理我还是建议用 Python 方式操作这也是多数气候数据处理流程中最省事的一环。2.2 用 Python GDAL 把温度读成 ndarrayGDAL 的 Grid 驱动可以自动识别这种目录式数据。假设你拥有一个Tem1991目录里面放着上述文件用下面的代码就能读到矩阵from osgeo import gdal import numpy as np path Tem1991 ds gdal.Open(path, gdal.GA_ReadOnly) if ds is None: raise IOError(无法打开栅格数据检查目录结构是否完整) print(行列数:, ds.RasterXSize, ds.RasterYSize) print(波段数:, ds.RasterCount) print(投影:, ds.GetProjection()) print(仿射变换参数:, ds.GetGeoTransform()) band ds.GetRasterBand(1) arr band.ReadAsArray() print(原始数组dtype:, arr.dtype, 最小值:, arr.min(), 最大值:, arr.max())这段代码的核心是gdal.Open注意传入的路径是Tem1991这个目录而不是arc0000.dat。如果路径正确GDAL 会自动识别 Grid 驱动。GetGeoTransform()返回六个值的元组分别是左上角 X 坐标、像元宽度、旋转项通常为 0、左上角 Y 坐标、旋转项、像元高度通常为负值。有了这组参数之后就可以把像元行列号转换成经纬度或投影坐标。很多人在这一步会直接读取.dat文件导致失败或者得到一维二进制流所以一定要从目录层面打开。2.3 数值还原与 NoData 处理数据说明里明确指出所有值都扩大了 10 倍除以 10 才是实际温度。因此读取后必须做一次换算# 将原始值转为浮点型并统一除以10得到摄氏度 temp_celsius arr.astype(np.float32) / 10.0 print(温度范围℃:, np.nanmin(temp_celsius), np.nanmax(temp_celsius))这里有一个容易被忽略的问题原始数组如果是整型除以 10 后得到的是小数所以要先astype(np.float32)再做除法否则整数除以整数会丢失小数位。另一个关键是 NoData 值。有些栅格在无数据区域使用-9999或0表示缺测如果不处理直接除以 10 后会得到无效的负值或 0干扰后续统计。建议在读取波段后立即检查 NoData 值nodata band.GetNoDataValue() print(NoData 原始值:, nodata) if nodata is not None: temp_celsius[temp_celsius nodata / 10.0] np.nan注意这里用nodata / 10.0与换算后的数组比较因为原始数组还没换算时NoData 值是原始编码。如果在换算前标记就直接比较原始数组然后赋 NaN 再换算。养成这个习惯后续做区域均值和趋势分析时才不会出现异常低值。降水数据同理只是单位从摄氏度变成毫米同样也要注意放大 10 倍的问题。3. 裁剪、重投影与像元归一把原始栅格变成分析可用数据3.1 为什么需要裁剪和重投影成渝城市群是一个行政区划概念而这份栅格数据很可能覆盖的是整个西南地区甚至全国。直接使用全区域计算不仅效率低而且边界外像元会干扰城市群尺度的统计。裁剪势在必行。另一个问题是投影。从数据产出来看原始数据很可能采用了 Albers 等积投影或兰伯特投影以便保证面积不变。但如果你之后要跟站点数据、气象插值数据合并或者用seaborn、geopandas画图通常需要把投影统一到 WGS84 地理坐标系EPSG:4326上这样经纬度才是直观的坐标。需要注意重投影会改变像元值吗如果使用最近邻采样像元值不会插值改变但面积和像元大小会变化如果使用双线性插值则会对温度值做平滑这取决于你的用途。气候栅格分析一般都建议用最近邻或者先完成数值统计再做投影转换。3.2 用 Rasterio 按成渝城市群边界裁剪Rasterio 是 GDAL 的 Python 封装处理栅格裁剪比直接写 GDAL 更简洁。先备好一份成渝城市群的矢量边界比如chengyu_boundary.shp然后用rasterio.mask.mask完成裁剪import geopandas as gpd import rasterio from rasterio.mask import mask as rio_mask boundary gpd.read_file(chengyu_boundary.shp) # 确保边界文件与栅格的CRS一致不一致时先转换 boundary boundary.to_crs(EPSG:4326) with rasterio.open(Tem1991) as src: # 读取原始投影若边界与栅格投影不同此处需要先把边界转到栅格投影 geoms boundary.geometry.values out_image, out_transform rio_mask(src, geoms, cropTrue, nodatasrc.nodata) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(Tem1991_chengyu.tif, w, **out_meta) as dst: dst.write(out_image)代码逻辑分三步第一步读取矢量边界第二步打开栅格并调用mask第三步把裁剪结果写为 GeoTIFF。cropTrue表示让输出范围严格贴合边界外接矩形nodata会继承原始栅格的 NoData 值。要注意的是mask函数要求边界的坐标系与栅格一致所以最好先把边界转换到栅格自身坐标系或者用rasterio.warp.transform_geom做几何转换。很多人在这一步报错就是因为投影不一致。另外如果你希望裁剪后的栅格保留原始投影而不是边界投影就不要提前把边界转成 EPSG:4326而是让 GDAL 自己处理。3.3 栅格统计区域均值怎么算才对裁剪完成后计算成渝城市群范围内的年均温均值最简单的办法是把像元值取平均但不要直接对裁剪后的数组做np.mean因为边界外的 NoData 会变成 0 或者无效值。正确姿势是使用np.nanmean并结合 NoData 掩膜。import numpy as np with rasterio.open(Tem1991_chengyu.tif) as src: band src.read(1, maskedTrue) # maskedTrue 会把NoData自动转为masked # band现在是一个masked array计算平均值时自动忽略无效值 region_mean band.mean() print(1991年成渝城市群平均气温由放大10倍数值还原:, float(region_mean) / 10.0, ℃)read(1, maskedTrue)返回一个numpy.ma.MaskedArray所有 NoData 都被隐藏在 mask 中mean()会跳过它们。如果不用 masked 模式就需要自己构建掩膜data src.read(1) nodata src.nodata valid data ! nodata mean_val data[valid].astype(np.float32).mean() / 10.0这两种方式都可以但maskedTrue在大规模计算时更直观。计算区域均值时还有一个容易被忽视的因素像元面积。如果栅格是等经纬度投影高纬度像元对应的面积比低纬度小此时计算区域平均温度应当做面积加权。但在成渝城市群这样一个南北跨度约 3 到 4 度的范围内面积差异很小做不等权校正的影响极其有限。如果你处理的是全国尺度数据建议用xarray的cos(latitude)加权。这里先按普通均值处理后续趋势分析也沿用这一口径确保不同年份之间统计方法一致。4. 1980—2015 年气温时空趋势斜率、显著性检验与可视化4.1 像元级线性趋势计算拿到 36 年年均温栅格后最想回答的问题通常是成渝城市群哪里在升温升温最剧烈的地方在哪这就需要对每个像元做时间序列回归。假设你已经把 1980 到 2015 年全部裁剪后的栅格读入一个三维数组形态为(36, rows, cols)那么逐像元的线性趋势可以用最小二乘拟合来求。为了效率不要用 Python 双重循环遍历所有像元而是把空间维展开成二维一次性计算years np.arange(1980, 2016, dtypenp.float64) # stack 形状:(36, rows, cols) # 转成 (36, rows*cols) n_years, rows, cols stack.shape flat stack.reshape(n_years, -1) # (36, N) # 中心化年份减少数值计算误差 years_c years - years.mean() denom (years_c ** 2).sum() # 计算每个像元的斜率 # 斜率 Σ( (year - year_mean) * (value - value_mean) ) / Σ(year - year_mean)^2 # 先计算每个像元在时间维上的方差 slope np.full(rows*cols, np.nan, dtypenp.float64) valid_mask np.all(~np.isnan(flat), axis0) if valid_mask.sum() 0: xn years_c.reshape(-1, 1) - years_c.mean()上面的写法有点绕更直接的方式是用向量化矩阵乘法# 标准化年份向量 x years - years.mean() x_std x / np.sqrt((x ** 2).sum()) # 对时间维中心化后的数据做回归 flat_centered flat - np.nanmean(flat, axis0, keepdimsTrue) # 斜率 sum(x_std * centered) / sqrt(sum(x^2)) , 由于x_std已经归一化 # 实际手动实现 y_demeaned flat - flat.mean(axis0, keepdimsTrue) # (36,N) slope_vec (x_std.reshape(-1,1) * y_demeaned).sum(axis0) / np.sqrt((x ** 2).sum())最终slope_vec就是每个像元的气温年倾率单位为℃/年。注意这里的数据必须已经是真温度也就是除以 10 后的数值。如果直接用原始编码值得到的斜率是 0.2 ℃/年 的 10 倍容易误判。计算完成后把slope_vecreshape 回二维就得到升温速率空间分布图。这个斜率就是经典的 Sen 斜率的一种简单线性估计适用没有异常值的年度序列。如果数据中有个别年份缺测np.isnan的掩膜会过滤掉但需要注意至少要有 30 年以上有效数据否则拟合结果稳健性差。4.2 区域平均温度序列的趋势检验像元级趋势能告诉你空间分布但政策或研究报告通常需要一条成渝城市群年平均气温区域序列并判断其变化是否统计显著。这里推荐 Mann-Kendall 趋势检验它不要求数据正态分布对异常值稳健。Python 中可用pymannkendall库实现pip install pymannkendall然后计算区域平均序列import pymannkendall as mk # 对每年stack取有效像元均值 annual_mean [] for i in range(stack.shape[0]): valid stack[i] annual_mean.append(np.nanmean(valid)) annual_mean np.array(annual_mean) # 单位℃ result mk.original_test(annual_mean, alpha0.05) print(MK趋势统计量:, result.Tau) print(p值:, result.p) print(趋势方向:, result.trend) print(斜率估计:, result.slope)pymannkendall的original_test返回一个结果对象其中trend会给出 increasing、decreasing 或 no trendslope则是 Sen 斜率单位是℃/年。相比普通最小二乘回归MK 检验不依赖线性假设能更保守地判断是否真的存在单调趋势。实测中成渝城市群年均温序列在 1980—2015 年间多呈现弱升温趋势但不同城市群子区域差异很大。如果你想分析四季或月平均气温趋势需要先把原始数据按季度合成再重复上述流程。4.3 趋势图与结果导出趋势结果最好用图来呈现。用matplotlib画两个图左图是像元趋势空间分布右图是区域平均时间序列与拟合线。先用basemap或cartopy绘制底图没有现成的城市群边界时直接使用裁剪后栅格的经纬度经纬度范围。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature lon np.linspace(transform[0], transform[0] transform[1]*cols, cols) lat np.linspace(transform[3] transform[5]*rows, transform[3], rows) fig plt.figure(figsize(12,5)) ax1 fig.add_subplot(1,2,1, projectionccrs.PlateCarree()) im ax1.imshow(slope_2d, extent[lon[0], lon[-1], lat[-1], lat[0]], cmapRdBu_r, transformccrs.PlateCarree()) ax1.add_feature(cfeature.COASTLINE) ax1.coastlines() plt.colorbar(im, axax1, shrink0.6, label℃/年) ax2 fig.add_subplot(1,2,2) ax2.scatter(years, annual_mean, s20, ck, alpha0.6) ax2.plot(years, slope* years intercept, r-, label线性拟合) ax2.set_xlabel(年份); ax2.set_ylabel(年均温℃) plt.show()绘制趋势图时要注意cartopy的imshow坐标参数extent需要按照左下角、右下角、右上角、左上角的经纬度顺序传入否则图像会上下颠倒。lat[-1]到lat[0]的写法是为了适应纬度递减的情况。趋势斜率图建议使用RdBu_r调色板红色代表升温蓝色代表降温符合气候报告的习惯。导出结果时可以直接用numpy.savetxt保存趋势矩阵或者用rasterio写成 GeoTIFF这样后续在 ArcGIS 中可以直接叠加查看。5. 批量转 NetCDF把 36 年温度数据装进一个文件5.1 用 xarray 合并逐年栅格逐年处理 36 个 GeoTIFF 非常麻烦而且每次读取都要重新打开文件。最理想的方式是把所有年份打包成一个多维 NetCDF维度为(time, latitude, longitude)之后用xarray一个文件搞定。首先要准备一个文件列表按年份排序然后逐个读取并构建坐标。import xarray as xr import numpy as np from osgeo import gdal file_list [fTem{year}_chengyu.tif for year in range(1980, 2016)] stack_list [] for year, path in zip(range(1980, 2016), file_list): ds gdal.Open(path) arr ds.ReadAsArray().astype(np.float32) / 10.0 gt ds.GetGeoTransform() # 计算每个像元的中心坐标注意gt[5]通常为负 lon gt[0] np.arange(ds.RasterXSize) * gt[1] gt[1] / 2 lat gt[3] np.arange(ds.RasterYSize) * gt[5] gt[5] / 2 stack_list.append(arr) ds None data np.array(stack_list) # (36, rows, cols)这里要注意gt[5]是负值表示纬度从北向南递减所以计算出的lat数组是递减的。xarray要求坐标单调递增因此需要翻转纬度轴data data[:, ::-1, :] # 翻转lat维度, 使其从小到大 lat lat[::-1] da xr.DataArray( data, dims(time, latitude, longitude), coords{time: list(range(1980, 2016)), latitude: lat, longitude: lon}, nametemperature, ) da.attrs[units] degC da.attrs[scale] original_values_divided_by_10xarray.DataArray在构建时如果坐标是浮点且不严格单调会自动报错。翻转后latitude变成从小到大的递增序列符合 CF 约定。如果你同时处理降水可以增加一个precipitation变量或者使用xr.Dataset包含两个变量。5.2 一次写出的 NetCDF 怎么用da.to_netcdf(chengyu_t2m_1980_2015.nc, enginenetcdf4)这行代码就会生成一个约几十兆的 NetCDF 文件包含全部 36 年的温度数据。之后再做任何分析都不需要再碰原始目录ds xr.open_dataset(chengyu_t2m_1980_2015.nc) temp ds[temperature] # (time, lat, lon) # 快速计算2000年之后的平均温度 recent_mean temp.sel(timeslice(2000, 2015)).mean(dimtime) # 快速提取某个点的温度序列 site_series temp.sel(longitude104.06, latitude30.67, methodnearest)xarray的切片、聚合、重采样能力比裸 NumPy 高效得多尤其是做多年平均、季节平均时直接用.groupby(time.year)或.resample(timeYS)即可。另外NetCDF 文件自带自描述信息在 Python 中读取不需要关心 NoDataxarray会将缺测值表示为NaN统计时会自动跳过。如果你要把数据交给 R 的ncdf4包处理或 Fortran 程序读取NetCDF 也是通用性最好的载体。这一步骤完成后原始数据的 10 倍放大问题已经彻底抹平后续所有分析的输入都是标准的物理量。本文还有配套的精品资源点击获取