企业官网建设流程全解析

做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. 工具选型：我为什么最后锁定了GDAL

Python里处理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:32650（WGS 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, dstSRS=dst_srs, format="GTiff", resamplingAlg="bilinear", creationOptions=["COMPRESS=LZW"] ) ds = None print("重投影完成，输出文件：", dst_path)

这段代码里几个参数我解释一下。dstSRS是目标坐标系，可以直接传EPSG代码、WKT字符串或者Proj4字符串。resamplingAlg指定重采样算法，这里我用的bilinear（双线性插值），适合连续型栅格数据，比如温度和DEM，可以让数值过渡更平滑。如果是类别型数据（土地利用类型、植被分类），必须用nearest neighbor（最近邻），否则会插出“四不像”的类别值。

creationOptions里的COMPRESS=LZW是给输出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)

如果你的目标是想让分