核心内容摘要
当我们聊2026展厅设计推荐时，会想到深圳北京河北设计团队
从 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(anaconda
然后选择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: #
获取文件的身份证信息元数据 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}) #
读取实际的像素数据 # 如果是卫星影像通常读取前三个波段红、绿、蓝 # 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(
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 矩阵分析鲁棒模式 ---) #
基础统计分析 - 使用 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} 米) #
计算地形起伏度 relief max_height - min_height print(f地形起伏度 (相对高度): {relief:.2f} 米) except Exception as e: print(f数据统计出错可能是数据全为空: {e}) #
空间筛选掩膜操作 # 注意直接比较 nan 1000 会得到 False通常不会报错但为了保险我们可以先处理一下 # 这里我们统计海拔大于 0 米且小于 9000 米排除异常值的像素 # 使用 符号连接两个条件 valid_elevation_mask (elevation_data
(elevation_data
high_altitude_mask (elevation_data
(elevation_data
# 计算符合条件的像素数量 pixel_count np.sum(high_altitude_mask) print(f海拔 1000米 的像素点数量: {pixel_count} 个) #
面积估算 pixel_area_sqm 30 * 30 total_high_area_km2 (pixel_count * pixel_area_sqm) / 1000000 print(f估算海拔 1000米 的区域面积: {total_high_area_km2:.2f} 平方公里) #
数据分布预览直方图数据的文本版 # 让我们看看数据大概集中在哪里使用 np.nanpercentile 查看分位数 # 这表示50% 的像素点海拔低于多少90% 的像素点海拔低于多少 p50 np.nanpercentile(elevation_data,
p90 np.nanpercentile(elevation_data,
print(f中位数海拔 (50%): {p50:.2f} 米) print(f90% 的区域海拔低于: {p90:.2f} 米) print(--- 分析结束 ---)运行代码后便可以获得如图所示的结果图7 分析结果掌握了对矩阵的运算你就拥有了分析数据的能力。
但光看枯燥的数字比如“平均海拔 500 米”往往不够直观。
接下来我们将使用 Matplotlib把 NumPy 矩阵重新变回可视化的图像绘制一张带有地理坐标的专业热力图。
使用 Matplotlib 绘制带地理坐标的热力图把数字变成颜色是可视化的重要形式。
在 Matplotlib 中这个过程叫做颜色映射。
我们可以利用代码制定一个规则数值越小低海拔涂成蓝色或绿色数值越大高海拔涂成棕色或白色。
Matplotlib 会自动根据我们提供的 DEM 矩阵中的数值大小给每一个像素点填上对应的颜色。
对于地形数据Python 还内置了一个专门的配色方案叫做 terrain地形它能自动模拟出从深海到平原、再到雪山的自然过渡色。
但仅仅画出颜色是不够的如果直接使用默认设置绘图会发现 X 轴和 Y 轴显示的数字是 0, 500,
..。
这代表的是像素的行列号而不是地理坐标。
对于专业的地理分析来说我们需要知道地图上某一点的具体经纬度。
这就涉及到了我们在之前提到的“Geo”信息。
为了给地图加上正确的坐标轴我们需要利用 Rasterio 读取到的边界信息。
这个边界信息告诉了我们这张图最左边是多少度最右边是多少度最上面和最下面又是多少度。
在绘图时我们将这些地理边界传递给 Matplotlib 的范围参数。
这样一来Matplotlib 就会明白这张图不是画在 0 到 2000 的像素坐标系里而是画在例如“东经
1
5 度到
1
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: #
读取 DEM 数据 (第1波段) data src.read(
#
处理异常值 # 将 -9999 等无效值设为 NaN这样绘图时这些区域就是透明的 data data.astype(float) data[data src.nodata] np.nan # 如果元数据没定义 nodata但你知道无效值是 -9999也可以手动写 # data[data -9999] np.nan #
获取地理边界 (关键步骤) # bounds 属性包含了 (左, 下, 右, 上) 的坐标边界 # 这决定了图片在坐标轴上的位置 spatial_extent (src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top) print(f数据地理范围: {spatial_extent}) # --- 开始绘图 --- # 设置画布大小单位是英寸dpi是分辨率 plt.figure(figsize(10,
, dpi
# 使用 imshow 绘制热力图 # cmapterrain: 使用内置的地形配色 (蓝-绿-棕-白) # extentspatial_extent: 将像素坐标映射到地理坐标 # vmin/vmax: 我们可以手动设置颜色映射的范围忽略极端的异常点 plt.imshow(data, cmapterrain, extentspatial_extent) # 添加颜色条 (Legend) # label 参数给颜色条加个说明 cbar plt.colorbar(labelElevation (meters)) # 添加标题和坐标轴标签 plt.title(Digital Elevation Model (DEM) Visualization, fontsize
plt.xlabel(Longitude (Degrees), fontsize
plt.ylabel(Latitude (Degrees), fontsize
# 添加网格线让坐标更清晰 plt.grid(True, linestyle--, alpha
0.
# 显示图像 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 #
创建基础网格 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**
#
模拟“红光波段” (Red Band) # 植物吸收红光所以森林区域(中心)红光值应该较低 # 我们设定中心低(
边缘高(
并加入随机噪点模拟真实环境 red_band 100 * distance_from_center np.random.randint(0, 20, (height, width)) # 限制数值范围在合理区间 red_band np.clip(red_band, 10,
#
模拟“近红外波段” (NIR Band) # 植物强烈反射近红外光所以森林区域(中心)值应该很高 # 我们设定中心高(
边缘低(
nir_band 200 * (1 - distance_from_center) np.random.randint(0, 20, (height, width)) nir_band np.clip(nir_band, 50,
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(divideignore, invalidignore) # 核心计算公式 ndvi (nir_float - red_float) / (nir_float red_float) print(fNDVI 计算完成。
) print(fNDVI 最大值: {np.nanmax(ndvi):.2f} (代表茂密植被)) print(fNDVI 最小值: {np.nanmin(ndvi):.2f} (代表裸土或水体)) # --- 步骤 3: 可视化对比 --- print(\n--- 步骤 3: 绘制对比图 ---) plt.figure(figsize(15,
, dpi
# 子图1: 红光波段 (看起来中心黑四周亮) plt.subplot(1, 3,
plt.imshow(red_band, cmapReds) plt.title(Simulated Red Band\n(Vegetation Absorbs Red)) plt.colorbar(labelPixel Value) # 子图2: 近红外波段 (看起来中心亮四周暗) plt.subplot(1, 3,
plt.imshow(nir_band, cmapGreys_r) # 使用灰度图亮色代表高值 plt.title(Simulated NIR Band\n(Vegetation Reflects NIR)) plt.colorbar(labelPixel Value) # 子图3: NDVI 结果 plt.subplot(1, 3,
# 使用 RdYlGn (红-黄-绿) 色带这是 NDVI 的标准配色 # vmin-
2, vmax
8 是 NDVI 的常见显示范围 plt.imshow(ndvi, cmapRdYlGn, vmin-
2, vmax
0.
plt.title(Calculated NDVI\n(Green Healthy Vegetation)) plt.colorbar(labelNDVI 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库安装完成然后我们再复制并运行以下代码import 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,
, dtypenp.uint
transform rasterio.transform.from_origin(
1
0,
3
0,
001,
0.
# 模拟经纬度坐标 # 写入临时文件 with rasterio.open( file_path, w, driverGTiff, heightdata.shape[0], widthdata.shape[1], count1, dtypedata.dtype, crsEPSG:4326, transformtransform ) as dst: dst.write(data,
print(f临时文件 {file_path} 已生成。
) # --- 正式开始裁剪流程 --- print(\n--- 开始执行空间裁剪 ---) #
定义裁剪区域模拟读取 Shapefile 的过程 # 我们创建一个矩形框 (minx, miny, maxx, maxy) # 假设我们要裁剪中间的一小块区域 # 原始范围大概是:
1
0 ~
1
5,
2
5 ~
3
0 minx, miny
1
1,
2
6 maxx, maxy
1
3,
2
8 # 使用 shapely 库构建一个几何对象 (GeoJSON 格式) bbox box(minx, miny, maxx, maxy) # rasterio 需要的格式是 GeoJSON 的 geometry 部分 shapes [bbox.__geo_interface__] print(f裁剪范围: 经度 {minx}~{maxx}, 纬度 {miny}~{maxy}) #
打开源文件并执行裁剪 with rasterio.open(file_path) as src: print(f源图像尺寸: {src.width} x {src.height}) # 检查坐标系 (实际项目中必须做的一步) # 如果 src.crs 和矢量数据的坐标系不一致这里需要进行投影转换 # 本例中我们手动生成的坐标是一致的所以跳过转换 # 执行裁剪 # cropTrue 表示把多余的空白边缘切掉只保留边界框内的数据 out_image, out_transform rasterio.mask.mask(src, shapes, cropTrue) # 更新元数据 (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]}) # ---
可视化对比 --- print(\n--- 绘制对比图 ---) plt.figure(figsize(12,
) # 显示原图 (为了方便显示我们需要重新读取原图) with rasterio.open(file_path) as src: original_data src.read(
plt.subplot(1, 2,
plt.imshow(original_data, cmapgray, 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
plt.title(Original Image with Clip Box) # 显示裁剪后的图 plt.subplot(1, 2,
# out_image 形状是 (1, h, w)绘图需要 (h, w) plt.imshow(out_image[0], cmapgray, extent(minx, maxx, miny, maxy)) plt.title(Clipped Result) plt.show() # ---
(可选) 保存裁剪结果 --- # with rasterio.open(clipped_result.tif, w, **out_meta) as dest: # dest.write(out_image) # print(裁剪结果已保存为 clipped_result.tif)运行代码后可以看到结果如下所示成功完成了裁剪图11 完成裁剪的效果到此我们便了解了如何从宏大的卫星影像中提取出重要的部分。
这不仅减少了数据量加快了后续的计算速度也让分析结果更加聚焦。
接下来我们开始探索地理数据处理的重投影问题。
重投影与重采样在实际工作中我们经常会遇到一个极其棘手的问题两份数据的坐标系不一样。
比如我们下载的卫星影像是 WGS84 (EPSG:
经纬度坐标但从测绘局拿到的行政区划图却是 CGCS2000 (EPSG:
投影坐标。
如果直接把它们叠在一起它们根本对不上。
这时候我们就需要进行重投影。
重投影把地球这个曲面上的图案换一种方式摊平到纸面上。
因为地球是圆的椭球体而地图是平的所以没有任何一种投影方式是完美的。
有的投影保面积如 Albers有的投影保角度如 Mercator。
当我们需要在不同投影之间转换数据时就是重投影。
而在重投影的过程中往往伴随着重采样。
例如你把一张正方形的网格图拉伸成菱形原本格子里的像素点位置就变了因此计算机需要决定新格子里的颜色该怎么填。
重采样有多种方法
最近邻法直接取最近那个格子的值。
适合分类数据如土地利用类型因为它不会产生新的数值不会把“森林”和“水体”平均成“沼泽”。
双线性插值取周围 4 个格子的平均值。
适合连续数据如卫星影像、气温会让画面更平滑。
三次卷积取周围 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: #
计算变换后的新参数 # 因为投影变了图像的宽、高、边界都会变 # 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}) #
创建目标文件并执行重投影 with rasterio.open(dst_filename, w, **kwargs) as dst: for i in range(1, src.count
: # 核心函数reproject # 它会把 src.read(i) 的数据填入 destination 数组中 reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, # 重采样算法选择 # 连续数据(如影像/高程)推荐用 Bilinear 或 Cubic # 分类数据(如土地利用)必须用 Nearest resamplingResampling.bilinear ) print(重投影完成) # ---
可视化对比 --- print(\n--- 绘制对比图 ---) plt.figure(figsize(12,
) # 读取原图 with rasterio.open(src_filename) as src: plt.subplot(1, 2,
plt.imshow(src.read(
, cmapgray) plt.title(fOriginal (EPSG:
\nUnits: Degrees) plt.axis(off) # 读取重投影后的图 with rasterio.open(dst_filename) as dst: plt.subplot(1, 2,
plt.imshow(dst.read(
, cmapgray) plt.title(fReprojected (EPSG:
\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。
这个过程中必须严格避免把小数存成整数否则
8 和
2 都会变成 0分析结果就会变成无意义状态。
复制并粘贴运行以下代码import rasterio import numpy as np # --- 准备数据模拟一个分析结果 --- print(正在准备待保存的数据...) ndvi_result np.random.uniform(-
5,
8, (500,
).astype(float
# --- 准备空间信息 --- # 在实际操作中这些信息通常来自原始影像 # src.transform 和 src.crs # 这里我们手动定义模拟一个位于北京附近的区域 # Transform 定义: (像素宽, 0, 左上角X, 0, 像素高(通常为负), 左上角Y) # 假设每个像素代表 30米 (约
0003 度) transform_info rasterio.transform.from_origin( west
1
3, # 左边界经度 north
4
0, # 上边界纬度 xsize
0003, # 像素宽度 (度) ysize
0003 # 像素高度 (度) ) crs_info rasterio.crs.CRS.from_epsg(
# WGS84 坐标系 output_file final_ndvi_analysis.tif print(f准备写入文件: {output_file}) print(f数据类型: {ndvi_result.dtype}) print(f坐标系: {crs_info}) # --- 正式写入流程 --- #
打开文件 (写入模式 w) # driverGTiff: 指定驱动为 GeoTIFF # height, width: 图像尺寸 # count: 波段数量 (NDVI 只有 1 个波段) # dtype: 数据类型 (非常重要必须与 ndvi_result.dtype 一致) # crs: 坐标系 # transform: 空间位置信息 # compresslzw: (可选) 使用 LZW 无损压缩可以减小文件体积 with rasterio.open( output_file, w, driverGTiff, heightndvi_result.shape[0], widthndvi_result.shape[1], count1, dtypendvi_result.dtype, crscrs_info, transformtransform_info, compresslzw ) as dst: #
写入数据 # 注意Rasterio 的 write 函数通常期望的形状是 (波段数, 高, 宽) # 或者我们可以指定写入到第几个波段 (indexes
# 如果 ndvi_result 是二维数组 (500,
直接写入会自动匹配到第1波段 dst.write(ndvi_result,
# (可选)
写入一些额外的描述信息 (Tags) dst.update_tags( CREATED_BYPython Rasterio Tutorial, DESCRIPTIONNDVI Analysis Result, DATE
) 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(
print(f读取到的数据均值: {np.mean(data_check):.4f}) print(验证通过该文件可以在 QGIS 或 ArcGIS 中正常打开。
)运行代码后观察结果可以发现我们已经成功保存了图13 完成保存的提示信息在用户文件夹中可以看到文件已经正常保存了图14 保存的文件至此我们便完整地走完了 Python 地理数据处理的全流程能够通过代码去挖掘地图背后隐藏信息了。