做GIS数据处理这段时间跟tiff栅格文件打交道的频率非常高。尤其是搞遥感影像、地形分析、气象插值这类活儿手里的数据来源杂得很有的坐标系是WGS84经纬度有的又是UTM投影有的分辨率是30米有的又是250米。数据之间互相不匹配做镶嵌、裁剪、波段运算之前必须先把它们统一到同一个坐标系、同一个像元尺寸下。这个统一的过程就是标题里说的重投影和重采样。我用Python跑这套流程已经两年多了从最开始拿命令行工具一个个敲到后来用GDAL、rasterio写批处理脚本折腾了不少弯路。这篇就把我用Python处理tiff栅格数据时关于重投影、重采样的完整思路和踩坑经验整理出来。如果你是做遥感、测绘、环境数据分析的或者刚接触栅格数据处理的这篇应该能帮你省下不少试错时间。1. 这俩操作到底在解决什么问题先搞清楚业务场景很多新手拿到tiff文件一上来就找代码跑重投影结果跑完发现结果不对又搞不清楚问题出在哪。我建议先停下来想清楚一件事你手里的数据到底为什么需要重投影或重采样1.1 重投影给数据换一套“定位语言”栅格数据的本质是一张由像元组成的二维矩阵每个像元记录一个数值比如高程、反射率、温度。这张矩阵要落到地球上必须有一套坐标系来告诉软件“每个像元对应地球表面的哪个位置”。这就是地理坐标系或投影坐标系的来历。问题在于不同数据源用的坐标系往往不一样。举个例子从USGS下载的Landsat影像默认是UTM投影而从某些气象站点下载的插值产品用的是WGS84经纬度。你把它们叠在一起看一个以米为单位、一个以度为经纬度单位位置完全对不上后续做任何叠加分析都是错的。重投影干的事就是把栅格数据从原来的坐标系“翻译”到目标坐标系。要注意的是这个翻译不只是坐标数值变了像元大小、图像形状都可能跟着变。因为同一块地面区域在不同的坐标系下表达出来的“网格铺法”不一样。这也是为什么重投影之后再检查数据你会看到分辨率、行数列数都变了。1.2 重采样让每个像元都有统一“尺寸”重采样解决的是另一个维度的问题像元尺寸不一致。你手头有一份30米分辨率的DEM还有一份500米分辨率的温度数据想算温度随高程的变化率就得先把两个栅格对到同一个网格上。要么把30米的降采样到500米要么把500米的升采样到30米但这两种做法代价完全不同。升采样听起来很美好——把小像元变细等于让数据“变清晰”了实际上这纯属错觉。重采样算法只能基于现有数值插值不可能无中生有地制造真实细节。你把500米的数据拉到30米只是把每个粗像元的值平滑地铺到更小的网格上伪增量没有新信息。降采样倒是合理的相当于把细网格的数据聚合到粗网格上会损失空间细节但保留了区域平均值这是做尺度上推时的标准操作。理解这一点很重要因为我见太多人一上来就把所有数据升采样到最高分辨率结果文件体积暴涨好几倍分析精度却没有实质提升。2. 工具选型我为什么最后锁定了GDALPython里处理tiff的工具其实不少rasterio、rioxarray、GDAL都能干这事。但它们之间的关系和各自定位新手容易搞混。我梳理一下自己的使用心得。2.1 Python生态里能干的活儿和各自的坑先说GDAL这是老牌地理数据抽象库几乎所有开源GIS软件的底层都在用它。它最大的优点是功能全、稳定重投影、重采样、格式转换、裁剪、镶嵌一个库全覆盖。缺点呢就是接口偏底层写法比较啰嗦而且安装容易出问题。rasterio是在GDAL之上封装的Python库API设计得非常Pythonic打开文件、读取数组、写入tiff都很顺手。它还没完全替代GDAL因为有些底层能力它没有完全暴露出来最终还是要回到GDAL的接口。rioxarray是基于xarray的栅格库处理多维栅格、带时间序列的遥感数据非常爽底层同样调用GDAL。如果你经常要处理NetCDF或多波段时序数据rioxarray很合适。如果只是做单张tiff的重投影和重采样拿它有点大材小用。我个人的选型经验是批处理脚本用GDAL的gdal.Warp和gdal.Translate配合os.listdir遍历文件速度快、代码直观做数据分析、要把栅格读成numpy数组来算的时候用rasterio打开因为它读出的数组和地理信息绑得很紧不容易搞错。2.2 GDAL的核心思路和装法GDAL最容易被忽略的点是它把读写操作统一抽象成了数据集Dataset和波段Band。一个tiff文件打开后是一个Dataset里面包含栅格的宽度、高度、地理变换参数、投影信息以及一个或多个Band。每个Band是二维数组存放不同波段的数据。重投影、重采样这些操作本质上是对Dataset层面的地理信息做修改同时重新组织像元数值。安装GDAL最稳妥的办法是用conda它能把相关依赖一并装好conda install -c conda-forge gdal如果非要用pip建议用预编译的wheel包否则源码编译会让你崩溃的。装完以后在Python里gdal.__version__验证一下能正常输出版本号就行了。3. 重投影实操从WGS84到UTM的完整流程直接上能跑的代码。我先说场景我从网上获取到一个WGS84坐标系下的tiff文件分辨率是0.01度约1公里但我的研究区在某个UTM分带内需要把数据转成UTM投影方便跟其他高分辨率数据叠加。目标投影我选择该区域对应的UTM分区这里我用EPSG:32650WGS 84 / UTM zone 50N。3.1 准备工作先看清数据原来的“身份证”动手之前先用GDAL把数据的投影、地理范围、像元大小这些基本信息摸清楚。这一步特别关键很多人连原数据是什么坐标系都没确认直接一顿操作最后结果自然没法看。from osgeo import gdal # 打开原始tiff src_path input_wgs84.tif ds gdal.Open(src_path) # 读取投影信息 src_prj ds.GetProjection() print(原始投影, src_prj) # 读取仿射变换参数 gt ds.GetGeoTransform() print(地理变换参数, gt) print(像元宽度, gt[1]) print(像元高度, gt[5]) # 读取栅格尺寸 print(列数, ds.RasterXSize, 行数, ds.RasterYSize) # 读取数据范围 min_x gt[0] max_y gt[3] max_x gt[0] gt[1] * ds.RasterXSize min_y gt[3] gt[5] * ds.RasterYSize print(数据范围, min_x, min_y, max_x, max_y)这里GetGeoTransform返回的六个参数依次是左上角X坐标、像元宽度、旋转项一般tiff里为0、左上角Y坐标、旋转项一般为0、像元高度注意是负值因为图像Y轴向下。这个数组是所有空间操作的基础务必理解。如果原数据的投影是“未知”的或者显示为空那说明tiff里压根没写坐标信息。这种情况没法直接重投影你得先从数据源确认它的坐标系再用gdal.Translate或SetProjection手动指定。这种问题我遇到很多次八成是数据在下载或格式转换时把坐标信息丢了。3.2 用gdal.Warp做投影转换GDAL做重投影的核心接口是gdal.Warp。它的含义很形象原本规则排列的像元网格要“扭曲”warp到新的坐标系下。gdal.Warp会根据源数据的像元值在新的网格上重新采样所以我们常说“重投影”其实天然就包含了“重采样”这一步。from osgeo import gdal src_path input_wgs84.tif dst_path output_utm50n.tif # 设置目标坐标系 dst_srs EPSG:32650 # 重投影 ds gdal.Warp( dst_path, src_path, dstSRSdst_srs, formatGTiff, resamplingAlgbilinear, creationOptions[COMPRESSLZW] ) ds None print(重投影完成输出文件, dst_path)这段代码里几个参数我解释一下。dstSRS是目标坐标系可以直接传EPSG代码、WKT字符串或者Proj4字符串。resamplingAlg指定重采样算法这里我用的bilinear双线性插值适合连续型栅格数据比如温度和DEM可以让数值过渡更平滑。如果是类别型数据土地利用类型、植被分类必须用nearest neighbor最近邻否则会插出“四不像”的类别值。creationOptions里的COMPRESSLZW是给输出tiff加无损压缩。栅格数据往往很占空间加压缩在文件大小上非常明显而且LZW是无损的不必担心数据失真。3.3 参数设置与验证重投影跑完不是生成文件就万事大吉了。我每次都会做两步验证一是看输出文件的投影信息是不是目标坐标系二是确认输出范围、分辨率是否符合预期。from osgeo import gdal result gdal.Open(dst_path) result_prj result.GetProjection() print(输出投影, result_prj) gt_new result.GetGeoTransform() print(新像元大小, gt_new[1], gt_new[5]) print(新列数, result.RasterXSize, 新行数, result.RasterYSize)如果你的目标是想让分
Python处理Tiff栅格数据：重投影与重采样实战教程
做GIS数据处理这段时间跟tiff栅格文件打交道的频率非常高。尤其是搞遥感影像、地形分析、气象插值这类活儿手里的数据来源杂得很有的坐标系是WGS84经纬度有的又是UTM投影有的分辨率是30米有的又是250米。数据之间互相不匹配做镶嵌、裁剪、波段运算之前必须先把它们统一到同一个坐标系、同一个像元尺寸下。这个统一的过程就是标题里说的重投影和重采样。我用Python跑这套流程已经两年多了从最开始拿命令行工具一个个敲到后来用GDAL、rasterio写批处理脚本折腾了不少弯路。这篇就把我用Python处理tiff栅格数据时关于重投影、重采样的完整思路和踩坑经验整理出来。如果你是做遥感、测绘、环境数据分析的或者刚接触栅格数据处理的这篇应该能帮你省下不少试错时间。1. 这俩操作到底在解决什么问题先搞清楚业务场景很多新手拿到tiff文件一上来就找代码跑重投影结果跑完发现结果不对又搞不清楚问题出在哪。我建议先停下来想清楚一件事你手里的数据到底为什么需要重投影或重采样1.1 重投影给数据换一套“定位语言”栅格数据的本质是一张由像元组成的二维矩阵每个像元记录一个数值比如高程、反射率、温度。这张矩阵要落到地球上必须有一套坐标系来告诉软件“每个像元对应地球表面的哪个位置”。这就是地理坐标系或投影坐标系的来历。问题在于不同数据源用的坐标系往往不一样。举个例子从USGS下载的Landsat影像默认是UTM投影而从某些气象站点下载的插值产品用的是WGS84经纬度。你把它们叠在一起看一个以米为单位、一个以度为经纬度单位位置完全对不上后续做任何叠加分析都是错的。重投影干的事就是把栅格数据从原来的坐标系“翻译”到目标坐标系。要注意的是这个翻译不只是坐标数值变了像元大小、图像形状都可能跟着变。因为同一块地面区域在不同的坐标系下表达出来的“网格铺法”不一样。这也是为什么重投影之后再检查数据你会看到分辨率、行数列数都变了。1.2 重采样让每个像元都有统一“尺寸”重采样解决的是另一个维度的问题像元尺寸不一致。你手头有一份30米分辨率的DEM还有一份500米分辨率的温度数据想算温度随高程的变化率就得先把两个栅格对到同一个网格上。要么把30米的降采样到500米要么把500米的升采样到30米但这两种做法代价完全不同。升采样听起来很美好——把小像元变细等于让数据“变清晰”了实际上这纯属错觉。重采样算法只能基于现有数值插值不可能无中生有地制造真实细节。你把500米的数据拉到30米只是把每个粗像元的值平滑地铺到更小的网格上伪增量没有新信息。降采样倒是合理的相当于把细网格的数据聚合到粗网格上会损失空间细节但保留了区域平均值这是做尺度上推时的标准操作。理解这一点很重要因为我见太多人一上来就把所有数据升采样到最高分辨率结果文件体积暴涨好几倍分析精度却没有实质提升。2. 工具选型我为什么最后锁定了GDALPython里处理tiff的工具其实不少rasterio、rioxarray、GDAL都能干这事。但它们之间的关系和各自定位新手容易搞混。我梳理一下自己的使用心得。2.1 Python生态里能干的活儿和各自的坑先说GDAL这是老牌地理数据抽象库几乎所有开源GIS软件的底层都在用它。它最大的优点是功能全、稳定重投影、重采样、格式转换、裁剪、镶嵌一个库全覆盖。缺点呢就是接口偏底层写法比较啰嗦而且安装容易出问题。rasterio是在GDAL之上封装的Python库API设计得非常Pythonic打开文件、读取数组、写入tiff都很顺手。它还没完全替代GDAL因为有些底层能力它没有完全暴露出来最终还是要回到GDAL的接口。rioxarray是基于xarray的栅格库处理多维栅格、带时间序列的遥感数据非常爽底层同样调用GDAL。如果你经常要处理NetCDF或多波段时序数据rioxarray很合适。如果只是做单张tiff的重投影和重采样拿它有点大材小用。我个人的选型经验是批处理脚本用GDAL的gdal.Warp和gdal.Translate配合os.listdir遍历文件速度快、代码直观做数据分析、要把栅格读成numpy数组来算的时候用rasterio打开因为它读出的数组和地理信息绑得很紧不容易搞错。2.2 GDAL的核心思路和装法GDAL最容易被忽略的点是它把读写操作统一抽象成了数据集Dataset和波段Band。一个tiff文件打开后是一个Dataset里面包含栅格的宽度、高度、地理变换参数、投影信息以及一个或多个Band。每个Band是二维数组存放不同波段的数据。重投影、重采样这些操作本质上是对Dataset层面的地理信息做修改同时重新组织像元数值。安装GDAL最稳妥的办法是用conda它能把相关依赖一并装好conda install -c conda-forge gdal如果非要用pip建议用预编译的wheel包否则源码编译会让你崩溃的。装完以后在Python里gdal.__version__验证一下能正常输出版本号就行了。3. 重投影实操从WGS84到UTM的完整流程直接上能跑的代码。我先说场景我从网上获取到一个WGS84坐标系下的tiff文件分辨率是0.01度约1公里但我的研究区在某个UTM分带内需要把数据转成UTM投影方便跟其他高分辨率数据叠加。目标投影我选择该区域对应的UTM分区这里我用EPSG:32650WGS 84 / UTM zone 50N。3.1 准备工作先看清数据原来的“身份证”动手之前先用GDAL把数据的投影、地理范围、像元大小这些基本信息摸清楚。这一步特别关键很多人连原数据是什么坐标系都没确认直接一顿操作最后结果自然没法看。from osgeo import gdal # 打开原始tiff src_path input_wgs84.tif ds gdal.Open(src_path) # 读取投影信息 src_prj ds.GetProjection() print(原始投影, src_prj) # 读取仿射变换参数 gt ds.GetGeoTransform() print(地理变换参数, gt) print(像元宽度, gt[1]) print(像元高度, gt[5]) # 读取栅格尺寸 print(列数, ds.RasterXSize, 行数, ds.RasterYSize) # 读取数据范围 min_x gt[0] max_y gt[3] max_x gt[0] gt[1] * ds.RasterXSize min_y gt[3] gt[5] * ds.RasterYSize print(数据范围, min_x, min_y, max_x, max_y)这里GetGeoTransform返回的六个参数依次是左上角X坐标、像元宽度、旋转项一般tiff里为0、左上角Y坐标、旋转项一般为0、像元高度注意是负值因为图像Y轴向下。这个数组是所有空间操作的基础务必理解。如果原数据的投影是“未知”的或者显示为空那说明tiff里压根没写坐标信息。这种情况没法直接重投影你得先从数据源确认它的坐标系再用gdal.Translate或SetProjection手动指定。这种问题我遇到很多次八成是数据在下载或格式转换时把坐标信息丢了。3.2 用gdal.Warp做投影转换GDAL做重投影的核心接口是gdal.Warp。它的含义很形象原本规则排列的像元网格要“扭曲”warp到新的坐标系下。gdal.Warp会根据源数据的像元值在新的网格上重新采样所以我们常说“重投影”其实天然就包含了“重采样”这一步。from osgeo import gdal src_path input_wgs84.tif dst_path output_utm50n.tif # 设置目标坐标系 dst_srs EPSG:32650 # 重投影 ds gdal.Warp( dst_path, src_path, dstSRSdst_srs, formatGTiff, resamplingAlgbilinear, creationOptions[COMPRESSLZW] ) ds None print(重投影完成输出文件, dst_path)这段代码里几个参数我解释一下。dstSRS是目标坐标系可以直接传EPSG代码、WKT字符串或者Proj4字符串。resamplingAlg指定重采样算法这里我用的bilinear双线性插值适合连续型栅格数据比如温度和DEM可以让数值过渡更平滑。如果是类别型数据土地利用类型、植被分类必须用nearest neighbor最近邻否则会插出“四不像”的类别值。creationOptions里的COMPRESSLZW是给输出tiff加无损压缩。栅格数据往往很占空间加压缩在文件大小上非常明显而且LZW是无损的不必担心数据失真。3.3 参数设置与验证重投影跑完不是生成文件就万事大吉了。我每次都会做两步验证一是看输出文件的投影信息是不是目标坐标系二是确认输出范围、分辨率是否符合预期。from osgeo import gdal result gdal.Open(dst_path) result_prj result.GetProjection() print(输出投影, result_prj) gt_new result.GetGeoTransform() print(新像元大小, gt_new[1], gt_new[5]) print(新列数, result.RasterXSize, 新行数, result.RasterYSize)如果你的目标是想让分
