rasterio实用教程（4）——坐标系转换

背景

在测绘领域中提到的坐标系有两种，一种是地理坐标系，用经纬度高程来表达；另一种是投影坐标系，即经过投影变换后的平面坐标系，通常是xy表达。

坐标系转换有四种情况：

1. 地理坐标系转投影坐标系（也叫大地坐标正算）
2. 投影坐标系转地理坐标系（也叫大地坐标反算）
3. 一种地理坐标系转另一种地理坐标系
4. 一种投影坐标系转另一种投影坐标系

本文的坐标系转换理论上涵盖以上四种情况。

实战

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject
from rasterio import crs

dataset = rasterio.open('D:/A_2021寒假/城市群相关/Data/ESACCI土地覆盖数据1992-2015/final/mask_proj_con_mask_ESACCI-LC-L4-LCCS-Map-300m-P1Y-1992-v2.0.7.tif')
src_img = 'input.tif' # 输入图像
dst_img = 'output.tif' # 输出图像
dst_crs = crs.CRS.from_proj4('+proj=lcc +lat_1=15 +lat_2=65 +lat_0=30 +lon_0=95 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs') # 输出图像坐标系

with rasterio.open(src_img) as src_ds:
	# 计算在新空间参考系下的仿射变换参数，图像尺寸
	dst_transform, dst_width, dst_height = calculate_default_transform(
		src_ds.crs, # 输入坐标系
		dst_crs, # 输出坐标系
		src_ds.width, # 输入图像宽
		src_ds.height, # 输入图像高
		*src_ds.bounds) # 输入数据源的图像范围

# 更新数据集的元数据信息
profile = src_ds.meta.copy()
profile.update({
	'crs': dst_crs,
	'transform': dst_transform,
	'width': dst_width,
	'height': dst_height
	})

# 重投影并写入数据
with rasterio.open(dst_img, 'w', **profile) as dst_ds:
	for i in range(1, src_ds.count + 1): # 遍历每个图层，通常只需要第一层即可
		src_array = src_ds.read(i)
		dst_array = np.empty((dst_height, dst_width), dtype=profile['dtype']) # 初始化输出图像数据
		
		# 重投影
		reproject(
			# 源文件参数
			source=src_array,
			src_crs=src_ds.crs,
			src_transform=src_ds.transform,
			# 目标文件参数
			destination=dst_array,
			dst_transform=dst_transform,
			dst_crs=dst_crs,
			# 其它配置
			resampling=Resampling.average,
			num_threads=2)
		# 写入图像
		dst_ds.write(dst_array, i)

如果你看过我上一篇博文重采样，你会发现代码基本一致。是的，实际上重采样和坐标系转换有个共同点，就是需要计算目标图像在新的坐标系下的仿射变换参数和空间范围。这一步都是通过calculate_default_transform函数来实现的。

坐标系转换的核心在于要写对输出图像的坐标系，即代码中的dst_crs，我这里采用的是crs.CRS.from_proj4(…)，实际上还有一些方法，例如from_wkt、from_epsg等，具体参考官网

如果你想获取正确的proj4 string，这里提供一个正规方法：

1. 打开spatial reference官网
2. 搜索你需要的坐标系，例如wgs84
3. 点击进入参考页面，如图
   
4. 点击Proj4，结果就出来了
   
5. 之后在from_proj4(…)函数中黏贴该字符串即可