Python 地理数据处理——GeoTIFF 读取与分析实战

Ne0inhk

23 Mar 2026 — 36 min read

一、从 TIF 到 GeoTIFF——当图像拥有了坐标

假设你用手机拍摄了一张风景照，或者在电脑上保存了一张高清的地图图片。在计算机眼里，这些普通的图片（通常是 JPG、PNG 或普通的 TIF 格式）其实只是由无数个像素排列而成的矩阵，计算机知道这张图有 1000 行、1000 列，知道第 5 行第 5 列的像素是绿色的，但它完全不知道这片绿色代表的是森林的一角，还是楼下的草坪。对于普通的图像处理软件来说，图片只是图片，它没有现实世界的概念，既不知道这张图在地球上的具体位置，也不知道图上的 1 厘米代表实际距离的 1 米还是 100 公里。

这就是普通 TIF 格式与 GeoTIFF 的核心区别所在。TIF本是一种为了存储高质量图像而设计的通用文件格式，常用于印刷和摄影。而 GeoTIFF，顾名思义，就是给这个普通的 TIF 文件穿上了一件“地理外衣”。它并没有改变图像原本的画面内容，而是在文件的内部“隐藏”了一系列特殊的标签信息。你可以把这些信息想象成贴在照片背面的一张详细说明书，上面写着：“这张照片的左上角对应地球上的北纬 30 度、东经 120 度”，“图片中的每一个像素点代表地面上 30 米 × 30 米的范围”，以及“这张图使用的是 WGS84 坐标系统”。

当 TIF 变成了 GeoTIFF，图像就拥有了空间意识。这种变化对于数据分析至关重要。因为有了坐标信息，我们不再是单纯地看图，而是可以进行运算。比如可以将一张 GeoTIFF 格式的卫星影像与另一张记录道路网络的矢量地图叠加在一起，只要它们拥有正确的地理坐标，无论这两份数据来源何处，它们都能在屏幕上严丝合缝地对齐。再比如，因为 GeoTIFF 告诉了计算机每一个像素代表的实际面积，就能直接通过数像素的方法，计算出某个湖泊的水域面积，或者某座城市的绿化覆盖率。

接下来，我们将不再使用普通的图片查看器，而是利用 Python 语言中的专业工具，去读取这些隐藏在图像背后的地理坐标信息，从而完成量化分析工作。

二、利用 Rasterio 加载卫星影像与 DEM 数据

在 Python 中，我们不能像双击图片那样直接查看内容，而是需要建立一个读取数据的通道。Rasterio 的工作流程非常符合这种逻辑：它使用一种被称为上下文管理器的语法，像是向图书馆借书——你借出书（打开文件），阅读内容，然后离开时系统会自动帮你把书还回去（关闭文件），防止文件被占用或损坏。

在使用 Rasterio 加载数据时，我们通常关注两部分内容：元数据和像素数据。当我们通过代码建立连接后，首先读取的不是巨大的图像本身，而是文件的元数据，包括了图像的宽度、高度、波段数量（比如卫星图通常有红绿蓝三个波段，而 DEM 只有一个高度波段），以及最重要的坐标系信息。确认这些信息无误后，我们才使用 read() 指令将真正的图像内容加载到内存中，变成计算机可以计算的数字矩阵。

打开开始菜单中的Anaconda(anaconda3)，然后选择Jupyter notebook：

图1 打开Jupyter Notebook

等待加载完毕，自动弹出网页。选择File——New——Console,点击Select，打开Python命令界面：

图2 打开Console

图3 点击Select

通过这种方式，我们便打开了 python 环境，做好了运行代码的准备。接下来，我们需要准备数据集。以示例数据集为例，将附件的数据集名字重命名为satellite_image.csv，然后将其粘贴到C盘的用户（User）文件夹内：

图4 粘贴到用户文件夹内

然后，我们需要安装 Rasterio 这个库，为导入做准备。复制并运行以下命令：

!pip install rasterio

出现以下提示后，证明安装已经完成：

图5 Rasterio 安装完成

最后，我们复制并粘贴以下代码，按Shift+回车运行：

import rasterio # 导入数据库文件 file_path = 'satellite_image.tif' print("正在尝试打开文件...") # 使用 with 语句打开文件，src 是我们对这个文件的称呼（Source的缩写） with rasterio.open(file_path) as src:          # 1. 获取文件的身份证信息（元数据）     print(f"文件名称: {src.name}")     print(f"图像尺寸: 宽 {src.width} x 高 {src.height} 像素")     print(f"波段数量: {src.count} 个")          # 查看地理坐标系（CRS），这告诉我们要把图贴在地球的哪个位置     print(f"坐标系统: {src.crs}")          # 查看仿射变换参数（Transform），它定义了像素与地理坐标的数学关系     print(f"空间变换参数:\n{src.transform}")     # 2. 读取实际的像素数据     # 如果是卫星影像，通常读取前三个波段（红、绿、蓝）     # read() 方法会将图像转换为数字矩阵     # 注意：Rasterio 的波段索引是从 1 开始的，而不是 0     if src.count >= 3:         # 读取全部波段         image_data = src.read()         print(f"已加载卫星影像数据，数据形状为: {image_data.shape}")         # 输出格式通常是 (波段数, 行数, 列数)     else:         # 如果是 DEM 数据（通常只有 1 个波段）         # 我们只读取第 1 个波段         elevation_data = src.read(1)         print(f"已加载 DEM 高程数据，数据形状为: {elevation_data.shape}")         # 输出格式通常是 (行数, 列数)     # 当代码运行出这个缩进块后，文件会自动关闭      print("文件读取完成，连接已关闭。")

如图所示，我们可以看到文件完成了读取：

图6 完成了文件的读取

三、基于 NumPy 的栅格矩阵运算

当我们把 GeoTIFF 加载进来后，得到的 elevation_data 变量，在专业术语中被称为 NumPy 数组。我们可以把它想象成一个超级强大的 Excel 表格，只不过这个表格没有行号和列号的显示，只有纯粹的数据。在之前的 DEM 数据中，表格里的每一个格子填写的都是海拔高度。利用 NumPy，我们可以瞬间对这几百万个数据进行数学运算，而不需要人工去一个个查看。

首先是统计分析。如果我们想知道这片区域的最高点有多高，或者最低点在哪里，如果用肉眼去图片上找颜色最亮或最暗的点，既费时又不准确。但对于 NumPy 来说，这只是毫秒级的计算。我们可以直接向这个矩阵提问：“你的最大值是多少？”、“平均海拔是多少？”。这种能力让我们能迅速掌握研究区域的地形概况，是进行任何地理分析的第一步。

其次是空间掩膜与筛选。如果你想研究“海拔 1000 米以上的区域”，在传统的纸质地图上，你需要拿着笔沿着等高线一点点描画。但在 NumPy 中，我们可以进行“条件筛选”。我们可以写一行代码，告诉计算机：“把所有小于 1000 的数字都遮住（或者标记为无效），只保留大于 1000 的数字”。这样，我们就瞬间提取出了高海拔区域。这种技术在地理分析中应用极广，比如筛选出适合种植某种作物的海拔区间，或者计算洪水淹没线以下的受灾区域。

最后我们还可以进行矩阵间的代数运算，虽然目前的 DEM 只有一个波段，但如果我们有两张不同年份的 DEM 数据（比如一张是 2010 年的，一张是 2020 年的），我们可以直接用“2020年的矩阵”减去“2010年的矩阵”。NumPy 会自动将对应位置的像素值相减。得到的新矩阵中，正数代表地形长高了（可能是填海造陆或火山喷发），负数代表地形降低了（可能是水土流失或挖掘）。这就是所谓的“变化检测”，其核心本质就是简单的矩阵减法。

复制并粘贴以下代码运行：

import numpy as np print("--- 开始 NumPy 矩阵分析（鲁棒模式） ---") # 1. 基础统计分析 - 使用 nan-safe 函数 # np.nanmax 会自动忽略矩阵中的 nan 值 try:     max_height = np.nanmax(elevation_data)     min_height = np.nanmin(elevation_data)     mean_height = np.nanmean(elevation_data)     print(f"区域最高海拔: {max_height:.2f} 米")     print(f"区域最低海拔: {min_height:.2f} 米")     print(f"区域平均海拔: {mean_height:.2f} 米")     # 2. 计算地形起伏度     relief = max_height - min_height     print(f"地形起伏度 (相对高度): {relief:.2f} 米") except Exception as e:     print(f"数据统计出错，可能是数据全为空: {e}") # 3. 空间筛选（掩膜操作） # 注意：直接比较 nan > 1000 会得到 False，通常不会报错，但为了保险我们可以先处理一下 # 这里我们统计海拔大于 0 米且小于 9000 米（排除异常值）的像素 # 使用 & 符号连接两个条件 valid_elevation_mask = (elevation_data > 0) & (elevation_data < 9000) high_altitude_mask = (elevation_data > 1000) & (elevation_data < 9000) # 计算符合条件的像素数量 pixel_count = np.sum(high_altitude_mask) print(f"海拔 > 1000米 的像素点数量: {pixel_count} 个") # 4. 面积估算 pixel_area_sqm = 30 * 30 total_high_area_km2 = (pixel_count * pixel_area_sqm) / 1000000 print(f"估算海拔 > 1000米 的区域面积: {total_high_area_km2:.2f} 平方公里") # 5. 数据分布预览（直方图数据的文本版） # 让我们看看数据大概集中在哪里，使用 np.nanpercentile 查看分位数 # 这表示：50% 的像素点海拔低于多少，90% 的像素点海拔低于多少 p50 = np.nanpercentile(elevation_data, 50) p90 = np.nanpercentile(elevation_data, 90) print(f"中位数海拔 (50%): {p50:.2f} 米") print(f"90% 的区域海拔低于: {p90:.2f} 米") print("--- 分析结束 ---")

运行代码后，便可以获得如图所示的结果：

图7 分析结果

掌握了对矩阵的运算，你就拥有了分析数据的能力。但光看枯燥的数字（比如“平均海拔 500 米”）往往不够直观。接下来我们将使用 Matplotlib，把 NumPy 矩阵重新变回可视化的图像，绘制一张带有地理坐标的专业热力图。

四、使用 Matplotlib 绘制带地理坐标的热力图

把数字变成颜色，是可视化的重要形式。在 Matplotlib 中，这个过程叫做颜色映射。我们可以利用代码制定一个规则：数值越小（低海拔），涂成蓝色或绿色；数值越大（高海拔），涂成棕色或白色。Matplotlib 会自动根据我们提供的 DEM 矩阵中的数值大小，给每一个像素点填上对应的颜色。对于地形数据，Python 还内置了一个专门的配色方案叫做 terrain（地形），它能自动模拟出从深海到平原、再到雪山的自然过渡色。

但仅仅画出颜色是不够的，如果直接使用默认设置绘图，会发现 X 轴和 Y 轴显示的数字是 0, 500, 1000...。这代表的是像素的行列号，而不是地理坐标。对于专业的地理分析来说，我们需要知道地图上某一点的具体经纬度。这就涉及到了我们在之前提到的“Geo”信息。

为了给地图加上正确的坐标轴，我们需要利用 Rasterio 读取到的边界信息。这个边界信息告诉了我们这张图最左边是多少度，最右边是多少度，最上面和最下面又是多少度。在绘图时，我们将这些地理边界传递给 Matplotlib 的范围参数。这样一来，Matplotlib 就会明白：这张图不是画在 0 到 2000 的像素坐标系里，而是画在例如“东经 113.5 度到 114.2 度”的地球表面上。

通过结合颜色映射与地理范围，我们就能生成一张既美观又精准的专业级热力图。这张图不仅能展示地形的高低起伏，还能让你把鼠标放在图上的任意位置，读出该点的经纬度坐标。

复制并粘贴运行以下代码：

import rasterio import matplotlib.pyplot as plt import numpy as np # --- 准备工作：我们需要重新读取文件以获取坐标范围信息 --- file_path = 'satellite_image.tif' print("正在读取数据并准备绘图...") with rasterio.open(file_path) as src:     # 1. 读取 DEM 数据 (第1波段)     data = src.read(1)          # 2. 处理异常值     # 将 -9999 等无效值设为 NaN，这样绘图时这些区域就是透明的     data = data.astype(float)     data[data == src.nodata] = np.nan     # 如果元数据没定义 nodata，但你知道无效值是 -9999，也可以手动写：     # data[data == -9999] = np.nan     # 3. 获取地理边界 (关键步骤！)     # bounds 属性包含了 (左, 下, 右, 上) 的坐标边界     # 这决定了图片在坐标轴上的位置     spatial_extent = (src.bounds.left, src.bounds.right,                       src.bounds.bottom, src.bounds.top)          print(f"数据地理范围: {spatial_extent}") # --- 开始绘图 --- # 设置画布大小，单位是英寸，dpi是分辨率 plt.figure(figsize=(10, 8), dpi=100) # 使用 imshow 绘制热力图 # cmap='terrain': 使用内置的地形配色 (蓝-绿-棕-白) # extent=spatial_extent: 将像素坐标映射到地理坐标 # vmin/vmax: 我们可以手动设置颜色映射的范围，忽略极端的异常点 plt.imshow(data, cmap='terrain', extent=spatial_extent) # 添加颜色条 (Legend) # label 参数给颜色条加个说明 cbar = plt.colorbar(label='Elevation (meters)') # 添加标题和坐标轴标签 plt.title("Digital Elevation Model (DEM) Visualization", fontsize=15) plt.xlabel("Longitude (Degrees)", fontsize=12) plt.ylabel("Latitude (Degrees)", fontsize=12) # 添加网格线，让坐标更清晰 plt.grid(True, linestyle='--', alpha=0.5) # 显示图像 print("正在渲染图像窗口...") plt.show()

运行代码后，我们就能看到生成的地形图，而且坐标轴上显示的不再是像素行数，而是真实的地理坐标。至此，我们便完成了一次专业的数据可视化。

图8 可视化分析结果

五、多波段计算——NDVI（归一化植被指数）实操

当我们看一张普通的彩色卫星照片时，我们看到的是红、绿、蓝三种颜色的混合。但在遥感科学中，卫星通常会拍摄更多的波段，其中一种叫做近红外波段（NIR）。植物的叶子非常独特，它们会吸收红光来进行光合作用，但会强烈反射近红外光。利用这个特性，可以总结出一个叫做 NDVI（归一化植被指数）的公式。它的原理很简单：用（近红外 - 红光）除以（近红外 + 红光）。结果越接近 1，说明植被越茂密。

我们可以利用 NumPy 模拟生成一份包含“红光”和“近红外”波段的虚拟数据，创建一个场景：图像中央是一片茂密的森林（近红外高，红光低），而周围是荒漠。

复制并粘贴运行以下代码：

import numpy as np import matplotlib.pyplot as plt print("--- 步骤 1: 模拟生成多光谱卫星数据 ---") # 设定图像大小 (500 x 500 像素) width, height = 500, 500 # 1. 创建基础网格 x = np.linspace(-1, 1, width) y = np.linspace(-1, 1, height) X, Y = np.meshgrid(x, y) # 计算每个点到中心的距离，用于生成圆形的"森林" distance_from_center = np.sqrt(X**2 + Y**2) # 2. 模拟“红光波段” (Red Band) # 植物吸收红光，所以森林区域(中心)红光值应该较低 # 我们设定：中心低(20)，边缘高(100)，并加入随机噪点模拟真实环境 red_band = 100 * distance_from_center + np.random.randint(0, 20, (height, width)) # 限制数值范围在合理区间 red_band = np.clip(red_band, 10, 150) # 3. 模拟“近红外波段” (NIR Band) # 植物强烈反射近红外光，所以森林区域(中心)值应该很高 # 我们设定：中心高(180)，边缘低(60) nir_band = 200 * (1 - distance_from_center) + np.random.randint(0, 20, (height, width)) nir_band = np.clip(nir_band, 50, 255) print("数据模拟完成。") print(f"红光波段形状: {red_band.shape}, 近红外波段形状: {nir_band.shape}") # --- 步骤 2: 计算 NDVI --- print("\n--- 步骤 2: 执行 NDVI 矩阵运算 ---") # 公式：NDVI = (NIR - Red) / (NIR + Red) # 注意：必须先转换为浮点数，否则整数除法会丢失精度 nir_float = nir_band.astype(float) red_float = red_band.astype(float) # 忽略除以0的警告（虽然模拟数据不太可能出现0） np.seterr(divide='ignore', invalid='ignore') # 核心计算公式 ndvi = (nir_float - red_float) / (nir_float + red_float) print(f"NDVI 计算完成。") print(f"NDVI 最大值: {np.nanmax(ndvi):.2f} (代表茂密植被)") print(f"NDVI 最小值: {np.nanmin(ndvi):.2f} (代表裸土或水体)") # --- 步骤 3: 可视化对比 --- print("\n--- 步骤 3: 绘制对比图 ---") plt.figure(figsize=(15, 5), dpi=100) # 子图1: 红光波段 (看起来中心黑，四周亮) plt.subplot(1, 3, 1) plt.imshow(red_band, cmap='Reds') plt.title("Simulated Red Band\n(Vegetation Absorbs Red)") plt.colorbar(label='Pixel Value') # 子图2: 近红外波段 (看起来中心亮，四周暗) plt.subplot(1, 3, 2) plt.imshow(nir_band, cmap='Greys_r') # 使用灰度图，亮色代表高值 plt.title("Simulated NIR Band\n(Vegetation Reflects NIR)") plt.colorbar(label='Pixel Value') # 子图3: NDVI 结果 plt.subplot(1, 3, 3) # 使用 RdYlGn (红-黄-绿) 色带，这是 NDVI 的标准配色 # vmin=-0.2, vmax=0.8 是 NDVI 的常见显示范围 plt.imshow(ndvi, cmap='RdYlGn', vmin=-0.2, vmax=0.8) plt.title("Calculated NDVI\n(Green = Healthy Vegetation)") plt.colorbar(label='NDVI Index') plt.tight_layout() plt.show()

运行代码生成的结果如下图所示：

图9 生成结果

可以发现：虽然红光图和近红外图只是黑白的深浅变化，但把它们相除之后，中间那团森林在 NDVI 图上瞬间变成了鲜艳的绿色，这就是遥感算法的效果。接下来，我们将学习如何对这些数据进行空间裁剪。

六、空间裁剪——利用 Shapefile 矢量边界裁剪 TIF

针对卫星影像，如果我们需要只剪出其中五角星形状的一小块，就需要利用空间裁剪。在计算机中，这个五角星形状通常存储为 .shp 或 GeoJSON 格式的矢量文件。这种文件不包含像素，只包含边界的坐标点。

在 Python 中，我们依然使用 Rasterio 库，但这次需要配合它的一个子模块 mask。首先需要读取矢量边界：使用 fiona 或 geopandas 库读取 Shapefile 文件，提取出多边形的几何形状（Geometry）。然后统一坐标系，确保“剪刀”和“桌布”必须在同一个坐标系下。如果一个是 WGS84（经纬度），另一个是 UTM（投影坐标），裁剪就会失败。所以必须先检查并转换坐标系。最后便可以执行裁剪，使用 rasterio.mask.mask 函数，传入影像和几何形状。它会返回两个东西：裁剪后的图像矩阵，以及一个新的仿射变换参数（因为裁剪后图像变小了，左上角的坐标也变了）。完成后将裁剪下来的数据和新的参数写入一个新的 GeoTIFF 文件中。

我们需要模拟创建一个矢量边界，即一个简单的矩形框，然后用它来裁剪我们在之前生成的模拟 NDVI 图像（或者任意加载的 TIF 数据）。在这之前，需要先运行以下命令，安装 shapely 库：

!pip install shapely

显示如下信息时，即代表库已经安装完成：

图10 shapely库安装完成

然后，我们再复制并运行以下代码：i

mport rasterio import rasterio.mask from shapely.geometry import box import numpy as np import matplotlib.pyplot as plt import json # --- 准备工作：我们需要一个待裁剪的影像 --- # 我们先创建一个临时的 GeoTIFF 文件 file_path = 'temp_source.tif' print("正在生成临时 GeoTIFF 文件用于演示...") # 创建一个 500x500 的随机数据 data = np.random.randint(0, 255, (500, 500), dtype=np.uint8) transform = rasterio.transform.from_origin(114.0, 30.0, 0.001, 0.001) # 模拟经纬度坐标 # 写入临时文件 with rasterio.open(     file_path, 'w', driver='GTiff',     height=data.shape[0], width=data.shape[1],     count=1, dtype=data.dtype,     crs='EPSG:4326', transform=transform ) as dst:     dst.write(data, 1) print(f"临时文件 {file_path} 已生成。") # --- 正式开始裁剪流程 --- print("\n--- 开始执行空间裁剪 ---") # 1. 定义裁剪区域（模拟读取 Shapefile 的过程） # 我们创建一个矩形框 (minx, miny, maxx, maxy) # 假设我们要裁剪中间的一小块区域 # 原始范围大概是: 114.0 ~ 114.5, 29.5 ~ 30.0 minx, miny = 114.1, 29.6 maxx, maxy = 114.3, 29.8 # 使用 shapely 库构建一个几何对象 (GeoJSON 格式) bbox = box(minx, miny, maxx, maxy) # rasterio 需要的格式是 GeoJSON 的 geometry 部分 shapes = [bbox.__geo_interface__] print(f"裁剪范围: 经度 {minx}~{maxx}, 纬度 {miny}~{maxy}") # 2. 打开源文件并执行裁剪 with rasterio.open(file_path) as src:     print(f"源图像尺寸: {src.width} x {src.height}")          # 检查坐标系 (实际项目中必须做的一步)     # 如果 src.crs 和矢量数据的坐标系不一致，这里需要进行投影转换     # 本例中我们手动生成的坐标是一致的，所以跳过转换          # 执行裁剪     # crop=True 表示把多余的空白边缘切掉，只保留边界框内的数据     out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)          # 更新元数据 (Metadata)     # 因为图像变小了，宽度、高度、变换参数都需要更新     out_meta = src.meta.copy()     out_meta.update({         "driver": "GTiff",         "height": out_image.shape[1],         "width": out_image.shape[2],         "transform": out_transform     }) print(f"裁剪完成！新图像尺寸: {out_image.shape[2]} x {out_image.shape[1]}") # --- 3. 可视化对比 --- print("\n--- 绘制对比图 ---") plt.figure(figsize=(12, 6)) # 显示原图 (为了方便显示，我们需要重新读取原图) with rasterio.open(file_path) as src:     original_data = src.read(1)     plt.subplot(1, 2, 1)     plt.imshow(original_data, cmap='gray', extent=(src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top))     # 画出裁剪框的红线     plt.plot([minx, maxx, maxx, minx, minx], [miny, miny, maxy, maxy, miny], 'r-', linewidth=2)     plt.title("Original Image with Clip Box") # 显示裁剪后的图 plt.subplot(1, 2, 2) # out_image 形状是 (1, h, w)，绘图需要 (h, w) plt.imshow(out_image[0], cmap='gray', extent=(minx, maxx, miny, maxy)) plt.title("Clipped Result") plt.show() # --- 4. (可选) 保存裁剪结果 --- # with rasterio.open("clipped_result.tif", "w", **out_meta) as dest: #     dest.write(out_image) # print("裁剪结果已保存为 clipped_result.tif")

运行代码后，可以看到结果如下所示，成功完成了裁剪：

图11 完成裁剪的效果

到此，我们便了解了如何从宏大的卫星影像中提取出重要的部分。这不仅减少了数据量，加快了后续的计算速度，也让分析结果更加聚焦。接下来，我们开始探索地理数据处理的重投影问题。

七、重投影与重采样

在实际工作中，我们经常会遇到一个极其棘手的问题：两份数据的坐标系不一样。比如，我们下载的卫星影像是 WGS84 (EPSG:4326) 经纬度坐标，但从测绘局拿到的行政区划图却是 CGCS2000 (EPSG:4527) 投影坐标。如果直接把它们叠在一起，它们根本对不上。这时候，我们就需要进行重投影。

重投影把地球这个曲面上的图案，换一种方式摊平到纸面上。因为地球是圆的（椭球体），而地图是平的，所以没有任何一种投影方式是完美的。有的投影保面积（如 Albers），有的投影保角度（如 Mercator）。当我们需要在不同投影之间转换数据时，就是重投影。而在重投影的过程中，往往伴随着重采样。例如，你把一张正方形的网格图拉伸成菱形，原本格子里的像素点位置就变了，因此计算机需要决定新格子里的颜色该怎么填。重采样有多种方法：

1.最近邻法：直接取最近那个格子的值。适合分类数据（如土地利用类型），因为它不会产生新的数值（不会把“森林”和“水体”平均成“沼泽”）。

2.双线性插值：取周围 4 个格子的平均值。适合连续数据（如卫星影像、气温），会让画面更平滑。

3.三次卷积：取周围 16 个格子的加权平均。计算慢，但图像质量最高。

在 Python 中，Rasterio 提供了非常方便的 reproject 函数来帮我们完成这项复杂的数学变换。为了演示效果，我们将把之前生成的模拟数据转换成一种常用的投影坐标系。你会发现，转换后图像的形状和像素尺寸都会发生变化。复制并运行以下代码：

import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling import numpy as np import matplotlib.pyplot as plt # --- 准备工作：使用生成的临时文件 --- src_filename = 'temp_source.tif' dst_filename = 'reprojected_result.tif' # 目标坐标系：EPSG:3857 (Web Mercator，单位是米) # 原始坐标系是 EPSG:4326 (WGS84，单位是度) dst_crs = 'EPSG:3857' print(f"正在准备将 {src_filename} 重投影到 {dst_crs} ...") try:     with rasterio.open(src_filename) as src:         # 1. 计算变换后的新参数         # 因为投影变了，图像的宽、高、边界都会变         # calculate_default_transform 会自动帮我们算好这些新参数         transform, width, height = calculate_default_transform(             src.crs, dst_crs, src.width, src.height, *src.bounds)         # 更新元数据         kwargs = src.meta.copy()         kwargs.update({             'crs': dst_crs,             'transform': transform,             'width': width,             'height': height         })         print(f"原始尺寸: {src.width} x {src.height}")         print(f"新尺寸: {width} x {height}")         # 2. 创建目标文件并执行重投影         with rasterio.open(dst_filename, 'w', **kwargs) as dst:             for i in range(1, src.count + 1):                 # 核心函数：reproject                 # 它会把 src.read(i) 的数据，填入 destination 数组中                 reproject(                     source=rasterio.band(src, i),                     destination=rasterio.band(dst, i),                     src_transform=src.transform,                     src_crs=src.crs,                     dst_transform=transform,                     dst_crs=dst_crs,                     # 重采样算法选择：                     # 连续数据(如影像/高程)推荐用 Bilinear 或 Cubic                     # 分类数据(如土地利用)必须用 Nearest                     resampling=Resampling.bilinear                 )                  print("重投影完成！")     # --- 3. 可视化对比 ---     print("\n--- 绘制对比图 ---")     plt.figure(figsize=(12, 6))     # 读取原图     with rasterio.open(src_filename) as src:         plt.subplot(1, 2, 1)         plt.imshow(src.read(1), cmap='gray')         plt.title(f"Original (EPSG:4326)\nUnits: Degrees")         plt.axis('off')     # 读取重投影后的图     with rasterio.open(dst_filename) as dst:         plt.subplot(1, 2, 2)         plt.imshow(dst.read(1), cmap='gray')         plt.title(f"Reprojected (EPSG:3857)\nUnits: Meters")         plt.axis('off')     plt.tight_layout()     plt.show() except FileNotFoundError:     print(f"错误：找不到文件 {src_filename}。请先运行上一章代码生成该文件。") except Exception as e: print(f"发生错误: {e}")

运行代码后，观察结果能够发现两张图看起来拉伸程度不一样。EPSG:4326 (原图)是因为经纬度格网在不同纬度其实不是正方形，直接显示时可能会觉得图像有点扁。EPSG:3857 (新图)是因为墨卡托投影在局部保持了形状，所以看起来可能更接近正方形。而且，新图的像素数量通常会发生变化，不再是严格的 500x500。掌握了重投影，无论数据来自哪个国家、哪个机构，我们都能把它们统一到一个坐标系下进行叠加分析。

图12 完成重投影与重采样的效果

接下来，我们学习如何将处理好、分析好、裁剪好、重投影好的数据，规范地写入并保存为一个标准的 GeoTIFF 文件，以便在其他软件中使用。

八、写入带有完整空间信息的 GeoTIFF 文件

在地理数据处理中，保存文件比读取文件要繁琐一些。普通的 JPG 图片只需要保存颜色，而 GeoTIFF 必须同时保存数据本身、坐标系、仿射变换参数。如果缺少了后两项，保存的结果就会变成一张普通的图片，在GIS软件中无法与真实地图对齐。在 Rasterio 中，写入文件的核心是 rasterio.open(..., 'w', ...) 模式。我们需要在打开文件时，就把这三要素作为参数传进去。这就像是在新建一个文档时，必须先设定好纸张大小、页边距和字体一样。另外，还有一个容易被忽视的细节：数据类型 (dtype)。如果结果是 NDVI（小数，范围 -1 到 1），就必须保存为 float32 或 float64。如果保存的是分类结果（1 代表森林，2 代表水体），为了节省空间，则应该保存为 uint8 或 int16。这个过程中，必须严格避免把小数存成整数，否则 0.8 和 0.2 都会变成 0，分析结果就会变成无意义状态。

复制并粘贴运行以下代码：

import rasterio import numpy as np # --- 准备数据：模拟一个分析结果 --- print("正在准备待保存的数据...") ndvi_result = np.random.uniform(-0.5, 0.8, (500, 500)).astype('float32') # --- 准备空间信息 --- # 在实际操作中，这些信息通常来自原始影像 # src.transform 和 src.crs # 这里我们手动定义，模拟一个位于北京附近的区域 # Transform 定义: (像素宽, 0, 左上角X, 0, 像素高(通常为负), 左上角Y) # 假设每个像素代表 30米 (约 0.0003 度) transform_info = rasterio.transform.from_origin(     west=116.3,      # 左边界经度     north=40.0,      # 上边界纬度     xsize=0.0003,    # 像素宽度 (度)     ysize=0.0003     # 像素高度 (度) ) crs_info = rasterio.crs.CRS.from_epsg(4326) # WGS84 坐标系 output_file = "final_ndvi_analysis.tif" print(f"准备写入文件: {output_file}") print(f"数据类型: {ndvi_result.dtype}") print(f"坐标系: {crs_info}") # --- 正式写入流程 --- # 1. 打开文件 (写入模式 'w') # driver='GTiff': 指定驱动为 GeoTIFF # height, width: 图像尺寸 # count: 波段数量 (NDVI 只有 1 个波段) # dtype: 数据类型 (非常重要！必须与 ndvi_result.dtype 一致) # crs: 坐标系 # transform: 空间位置信息 # compress='lzw': (可选) 使用 LZW 无损压缩，可以减小文件体积 with rasterio.open(     output_file,     'w',     driver='GTiff',     height=ndvi_result.shape[0],     width=ndvi_result.shape[1],     count=1,     dtype=ndvi_result.dtype,     crs=crs_info,     transform=transform_info,     compress='lzw' ) as dst:          # 2. 写入数据     # 注意：Rasterio 的 write 函数通常期望的形状是 (波段数, 高, 宽)     # 或者我们可以指定写入到第几个波段 (indexes=1)     # 如果 ndvi_result 是二维数组 (500, 500)，直接写入会自动匹配到第1波段     dst.write(ndvi_result, 1)          # (可选) 3. 写入一些额外的描述信息 (Tags)     dst.update_tags(         CREATED_BY="Python Rasterio Tutorial",         DESCRIPTION="NDVI Analysis Result",         DATE="2023-10-27"     ) print("文件写入成功！") # --- 验证环节 --- # 让我们重新读取刚刚保存的文件，看看信息是否正确 print("\n--- 验证保存的文件 ---") with rasterio.open(output_file) as src:     print(f"文件名: {src.name}")     print(f"尺寸: {src.width} x {src.height}")     print(f"坐标系: {src.crs}")     print(f"波段数: {src.count}")     print(f"左上角坐标: {src.bounds.left}, {src.bounds.top}")          # 读取一点数据看看     data_check = src.read(1)     print(f"读取到的数据均值: {np.mean(data_check):.4f}")     print("验证通过！该文件可以在 QGIS 或 ArcGIS 中正常打开。")

运行代码后，观察结果可以发现，我们已经成功保存了：