python对tif数据重投影

一、不同投影坐标系的区别

地理坐标系（Geographic Coordinate System, GCS）和投影坐标系（Projected Coordinate System, PCS）是两种常见的坐标系统，它们在表示地理信息时有着不同的方式。以下是它们的主要区别：

1. 定义方式：

示例：

示例：

2. 表示方式：

例如：WGS 84（纬度：40°N，经度：75°W）

例如：UTM 区域 33N（X: 500000, Y: 4649776）

3. 使用场景：

应用：GPS 定位、卫星数据、气象数据等。

应用：精确的地图制图、地籍调查、城市规划、地理分析等。

4. 转换关系：

5. 优缺点：

总结：

这两种坐标系各有优势，通常在实际应用中，根据不同需求使用不同的坐标系统。

二、EPSG投影代码

EPSG（European Petroleum Survey Group）是一个全球地理信息系统（GIS）领域中广泛使用的坐标参考系统（CRS）代码标准。每个坐标参考系统（例如，地理坐标系统、投影坐标系统等）都有一个唯一的 EPSG 代码。

常见的 EPSG 代码示例：

1. WGS 84 (全球定位系统标准)：

2. EPSG:4326：这是一种地理坐标系统，使用经度和纬度来表示位置，通常与 GPS 和许多地理信息系统（GIS）应用程序兼容。

3. UTM (Universal Transverse Mercator)：

4. EPSG:32633：这是 UTM 投影系统中的一个常用区域，适用于东经6°至12°的区域（例如，德国地区）。

5. EPSG:32733：这是 UTM 投影系统中的另一区域，适用于南半球的相同经度范围。

6. Web Mercator：

7. EPSG:3857：这是在线地图服务（如 Google Maps、OpenStreetMap 和 Bing Maps）常用的投影坐标系，广泛应用于互联网地图。

8. WGS 84 / Pseudo-Mercator：

9. EPSG:4326：这种坐标系将地球的经纬度数据表示为度，并且经常用在地图投影中。

常见 CRS 类型：

1. 地理坐标系统（Geographic Coordinate Systems）：使用经度和纬度来描述地球上的位置，单位通常是度（°）。

2. 例如：EPSG:4326（WGS84）、EPSG:4269（NAD83）等。

3. 投影坐标系统（Projected Coordinate Systems）：使用平面坐标来描述位置，通常使用米或英尺作为单位。

4. 例如：EPSG:3857（Web Mercator）、EPSG:32633（UTM区33N）等。

如何查询 EPSG 代码？

你可以通过以下几种方式查找特定的 EPSG 代码：

1. EPSG 官网：EPSG.io 提供了一个搜索工具，你可以通过输入坐标系统名称或相关信息来查找代码。

2. GIS 软件：在 QGIS 或 ArcGIS 等地理信息软件中，你可以查看和选择不同的 EPSG 代码。

在 Python 中使用 EPSG 代码：

如果你使用 rasterio 或 pyproj 等库来处理地理数据，可以通过 EPSG 代码来指定坐标系统。

例如，使用 pyproj 定义投影坐标系：

from pyproj import CRS

# 定义 WGS 84 坐标系统
crs = CRS.from_epsg(4326)  # WGS 84
print(crs)


import rasterio
from rasterio.warp import calculate_default_transform

# 打开原始影像
with rasterio.open('input.tif') as src:
    # 设置目标坐标系统为 EPSG:4326 (WGS 84)
    transform, width, height = calculate_default_transform(
        src.crs, 'EPSG:4326', src.width, src.height, *src.bounds)
    print(transform, width, height)

三、一个例子

这里我展示一个利用python实现对tif的坐标转换

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Nov 23 11:37:13 2024

@author: Chistrong wen
University of Stuttgart, Institute of Geodesy
"""


import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
import matplotlib.pyplot as plt
import numpy as np

# 打开原始影像
src_path = '/Users/chistrongwen/Downloads/SM_V3.0/SM_2016_10.tif'  # 输入影像的路径
dst_crs = 'EPSG:3857'  # 目标坐标系 (例如 WGS84)

with rasterio.open(src_path) as src:
    # 计算目标影像的坐标转换和新尺寸
    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
    })
    
    # 创建一个数组来存储重新投影后的影像
    reprojected_image = np.zeros((height, width), dtype=np.float32)
    
    # 将重新投影后的数据保存到新的文件中
    with rasterio.open('reprojected_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=dst_crs,
                resampling=Resampling.nearest)
    
    # 显示原始影像和重新投影后的影像
    with rasterio.open(src_path) as src:
        original_image = src.read(1)  # 读取第一个波段

    # 绘制图形：显示原始影像和重新投影后的影像
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
    
    ax1.imshow(original_image, cmap='viridis')
    ax1.set_title('Original image')
    
    # 显示重新投影后的影像
    with rasterio.open('reprojected_output.tif') as reprojected_src:
        reprojected_image = reprojected_src.read(1)
    
    ax2.imshow(reprojected_image, cmap='viridis')
    ax2.set_title('reprojection (EPSG:4326)')
    
    plt.show()

我们可以看到投影前后的区别

❤️欢迎点赞收藏❤️

作者：我是水怪的哥