如何为Numpy数组添加投影与地理变换后保存为TIFF文件

我从LiDAR表面导入了一个.tif文件到Python中：

Existing_Surface = 'Path_To_Tiff/Existing_Surface.tif'

使用GDAL检查了投影信息：

Original_Surface = gdal.Open(Existing_Surface)
projection_Original = Original_Surface.GetProjection()
geotransform_Original = Original_Surface.GetGeoTransform()

导入的影像具有正确的投影和地理变换。

随后我利用shp文件裁剪该影像，并将新影像转换为Numpy数组（示例数组如下）：

Existing_example_array = [[ 0, 0, 1, 0, 0, 0, 0], 
                          [ 0, 1, 1, 1, 0, 0, 0], 
                          [ 0, 1, 1, 1, 1, 0, 0], 
                          [ 1, 1, 1, 1, 1, 0, 0], 
                          [ 0, 1, 1, 1, 1, 0, 0], 
                          [ 0, 1, 1, 0, 0, 0, 0]]

对数组内的数据进行处理后，尝试导出该Numpy数组：

data_folder = Path("Path_to_Tiff/Surface.tif")
tifffile.imsave( data_folder, Existing_example_array)

tifffile成功导出了数组，但投影和地理变换信息丢失。请问有没有方法为数组恢复投影和地理变换并重新导出为TIFF文件？

tifffile仅负责处理像素数据，不支持地理空间元数据的读写，需要用GDAL或rasterio这类地理空间库来导出带投影和地理变换的TIFF。

方法一：使用GDAL导出 #

通过GDAL创建新的TIFF数据集，写入处理后的数组，再设置投影和地理变换：

from osgeo import gdal
from pathlib import Path

# 替换为你处理后的数组、投影、地理变换（若为裁剪影像，需用裁剪后的geotransform）
processed_array = Existing_example_array
proj = projection_Original
gt = geotransform_Original

output_path = Path("Path_to_Tiff/Surface_with_proj.tif")
rows, cols = processed_array.shape
# 根据原影像数据类型调整，比如float32用gdal.GDT_Float32
data_type = gdal.GDT_Byte

# 创建TIFF驱动
driver = gdal.GetDriverByName('GTiff')
# 创建数据集：路径、宽度、高度、波段数、数据类型
dataset = driver.Create(str(output_path), cols, rows, 1, data_type)

# 设置地理变换和投影
dataset.SetGeoTransform(gt)
dataset.SetProjection(proj)

# 写入数组数据
dataset.GetRasterBand(1).WriteArray(processed_array)

# 释放资源
dataset.FlushCache()
dataset = None

方法二：使用rasterio导出 #

rasterio的地理空间数据API更简洁直观：

import rasterio
from rasterio.transform import from_gdal
from pathlib import Path

processed_array = Existing_example_array
proj = projection_Original
gt = geotransform_Original

output_path = Path("Path_to_Tiff/Surface_with_proj.tif")

# 构建输出元数据
meta = {
    'driver': 'GTiff',
    'height': processed_array.shape[0],
    'width': processed_array.shape[1],
    'count': 1,
    'dtype': processed_array.dtype,
    'crs': proj,
    'transform': from_gdal(*gt),
}

# 写入数据
with rasterio.open(output_path, 'w', **meta) as dst:
    dst.write(processed_array, 1)

注意事项 #

• 如果是裁剪后的影像，不能直接使用原文件的geotransform，需要在裁剪过程中记录新的地理变换参数（比如用GDAL Warp或rasterio裁剪时获取输出的transform）。
• 确保导出数据类型与原影像一致，避免精度丢失或数据异常。

内容的提问来源于stack exchange，提问作者Patstro