python对tif数据重投影

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

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

1. 定义方式：

地理坐标系 (GCS)：
- 地理坐标系使用球面或椭球体上的经纬度（纬度和经度）来描述地球表面上的位置。
- 纬度（Latitude）表示地球赤道平面上某点与赤道的角度，范围为 -90° 到 +90°。
- 经度（Longitude）表示某点相对于本初子午线（零度经线）的角度，范围为 -180° 到 +180°。
- GCS 是基于一个地球模型（通常是一个椭球体或球体），例如 WGS 84 或 NAD83。
示例：
- WGS 84（EPSG:4326）
- NAD83（EPSG:4269）
投影坐标系 (PCS)：
- 投影坐标系则将地球的三维曲面转换为二维平面，使用平面坐标（通常是X和Y坐标）来描述位置。
- 投影坐标系通常使用某种数学投影方式（如墨卡托投影、UTM投影等）将球面上的坐标转换为平面坐标系上的点。
- 这类坐标系常常使用米、英尺或其他线性单位来表示距离。
示例：
- UTM（Universal Transverse Mercator，EPSG:32633）
- Web Mercator（EPSG:3857）

2. 表示方式：

地理坐标系：
- 坐标单位是角度（度，°），例如，经度和纬度。
- 坐标系是三维的，描述地球表面上的位置。
例如：WGS 84（纬度：40°N，经度：75°W）
投影坐标系：
- 坐标单位通常是米（m）或英尺（ft），并且是二维的坐标。
- 它通过投影方法将地球表面（球面）投影到平面上，通常会失真。
例如：UTM 区域 33N（X: 500000, Y: 4649776）

3. 使用场景：

地理坐标系：
- 地理坐标系适合用于全球范围内的应用，尤其是在需要精确表示位置时（例如 GPS 定位、全球范围的地图、气候模型等）。
- 由于地球是一个三维的球体，地理坐标系能够更好地反映地球表面的真实情况。
应用：GPS 定位、卫星数据、气象数据等。
投影坐标系：
- 投影坐标系适合用于局部区域的地图和精确测量，特别是在需要计算距离、面积等物理量时。
- 因为它是基于平面坐标，因此失真问题在大范围投影时会变得显著。为了减小失真，常常需要使用合适的投影方法和区域划分。
应用：精确的地图制图、地籍调查、城市规划、地理分析等。

4. 转换关系：

地理坐标系通常用于表示全球位置，而投影坐标系适用于局部区域，因为平面上的坐标系统更适合测量短距离。
从地理坐标系转换到投影坐标系（例如，从 WGS 84 转换到 UTM 坐标系）是常见的需求。转换时需要用到合适的投影算法，如墨卡托投影、横轴墨卡托投影等。

5. 优缺点：

地理坐标系的优点：
- 全球适用，不受局部区域影响。
- 适合全球性的数据表示。
- 缺点：
  - 不适合做精确测量。
  - 对于局部区域的数据处理和分析来说，坐标单位（角度）不够直观，转换成平面坐标系后可以避免失真。
投影坐标系的优点：
- 精确度较高，适用于局部区域的距离、面积和角度测量。
- 方便绘制和分析地图。
- 缺点：
  - 只适用于较小范围的区域，投影过程中存在失真问题。
  - 不适用于全球数据的表示。

总结：

地理坐标系：用于描述地球表面上的位置，坐标单位为经纬度（角度），适合全球范围内的数据。
投影坐标系：将地球表面的三维坐标转换为二维平面坐标，坐标单位通常为米或英尺，适合局部区域的数据处理和测量。

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

二、EPSG投影代码

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

常见的 EPSG 代码示例：

WGS 84 (全球定位系统标准)：
- EPSG:4326：这是一种地理坐标系统，使用经度和纬度来表示位置，通常与 GPS 和许多地理信息系统（GIS）应用程序兼容。
UTM (Universal Transverse Mercator)：
- EPSG:32633：这是 UTM 投影系统中的一个常用区域，适用于东经6°至12°的区域（例如，德国地区）。
- EPSG:32733：这是 UTM 投影系统中的另一区域，适用于南半球的相同经度范围。
Web Mercator：
- EPSG:3857：这是在线地图服务（如 Google Maps、OpenStreetMap 和 Bing Maps）常用的投影坐标系，广泛应用于互联网地图。
WGS 84 / Pseudo-Mercator：
- EPSG:4326：这种坐标系将地球的经纬度数据表示为度，并且经常用在地图投影中。

常见 CRS 类型：

地理坐标系统（Geographic Coordinate Systems）：使用经度和纬度来描述地球上的位置，单位通常是度（°）。
- 例如：EPSG:4326（WGS84）、EPSG:4269（NAD83）等。
投影坐标系统（Projected Coordinate Systems）：使用平面坐标来描述位置，通常使用米或英尺作为单位。
- 例如：EPSG:3857（Web Mercator）、EPSG:32633（UTM区33N）等。

如何查询 EPSG 代码？

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

EPSG 官网：EPSG.io 提供了一个搜索工具，你可以通过输入坐标系统名称或相关信息来查找代码。
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()

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