Python rasterio库进行投影转换

rasterio对遥感影像进行投影转换

在处理遥感数据时，常见的格式包括HDF、NC、DAT、TIF、GRIB等等，其中DAT、TIF数据常常涉及到投影坐标的问题，我们在进行不同源数据对比分析的时候，当两个数据源为不同投影类型时，无法直接进行分析，需要进行投影转换，使其空间上一一对应（像元尺度）。

1. 安装rasterio库

pip install rasterio
## 或者是用conda安装
conda install rasterio
## 安装失败的话  请检查一下是否打开了网络代理和VPN，需要在代理和VPN关闭状态下执行

2. 代码如下

from osgeo import gdal,ogr,osr
import scipy.io as sio
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
import rasterio
from rasterio.windows import from_bounds
from rasterio.transform import Affine
import numpy as np
## rasterio库好像是依赖gdal的，所以需要import一下osgeo
def reproject_to_wgs84(input_tif, output_tif):
    """将TIFF影像重投影为WGS-84 (EPSG:4326)"""
    # 打开源文件
    with rasterio.open(input_tif) as src:
        # 计算目标变换参数
        transform, width, height = calculate_default_transform(
            src.crs,               # 输入坐标系（自动读取）
            'EPSG:4326',           # 目标坐标系（WGS-84）
            src.width,             # 输入影像宽度
            src.height,            # 输入影像高度
            *src.bounds            # 输入影像边界 [left, bottom, right, top]
        )

        # 设置输出文件参数
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': 'EPSG:4326',    # 目标坐标系
            'transform': transform, # 新变换参数
            'width': width,         # 新宽度
            'height': height        # 新高度
        })

        # 执行重投影
        with rasterio.open(output_tif, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),      # 输入波段
                    destination=rasterio.band(dst, i), # 输出波段
                    src_transform=src.transform,       # 输入变换参数
                    src_crs=src.crs,                   # 输入坐标系
                    dst_transform=transform,           # 输出变换参数
                    dst_crs='EPSG:4326',               # 目标坐标系
                    resampling=Resampling.bilinear     # 重采样方法
                )
    print(f"重投影完成")
if __name__=="__main__":
	input_tif = '/user/original_img.tif'
	output_tif = '/user/out_img.tif'
	reproject_to_wgs84(input_tif, output_tif)

3.代码说明：
（1）输入数据为我们需要转换投影类型的数据“original_img.tif”
（2）输出数据为转投影后的结果，“out_img.tif”
（3）投影类型选择：EPSG:4326 代表wgs84坐标

4.转换结果如下：