【GDAL-Python】4-在Python中使用GDAL进行投影转换、影像重采样以及裁剪栅格数据

news2024/9/24 9:24:05

文章目录

1-介绍
- 1.1 主要内容
- 1.2 EPSG
2-代码实现
- 2.1 数据介绍
- 2.2 代码实现
- 2.3 结果及效果显示
3.参考资料

1-介绍

1.1 主要内容

（1）在本教程中，将介绍如何在 Python 中使用 gdalwarp 将栅格数据重新投影到不同的坐标参考系，影像重采样（改变影像分辨率）、以及根据矢量 shapefile文件将影像对应部分进行裁剪；
（2）技术路线：在python脚本中使用gdal.Warp，不同的功能对应不同的可选项参数

1）坐标参考系改变对应：gdal.Warp(dstSRS = “EPSG:4326”)目标参考系，采用EPSG方式。
2）影像分辨率gdal.Warp(“demRes.tif”, ds, xRes = 150, yRes = 150,resampleAlg），xRes水平采样分辨率，yRes垂直方向采样分辨率，resampleAlg对应采样的方法）；
3）gdal.Warp(“demClip.tif”, ds, cutlineDSName = “star.shp”,cropToCutline = True, dstNodata = np.nan)，对应参数cutlineDSNane支持矢量shp文件；

（3）视频地址：B站对应教程-4-在Python中使用GDAL进行投影转换、影像重采样以及裁剪栅格数据

1.2 EPSG

（1）EPSG代码是一种用于标准化地理空间参考系统（包括坐标系、投影、地理坐标系等）的唯一编码。这些编码帮助用户在不同的地理信息系统（GIS）软件和数据集中准确地识别和使用特定的地理空间参考系统。每个EPSG代码都对应着一组特定的参数，如椭球体、单位、地理坐标系或投影方式，这使得它们在全球范围内具有普遍的应用价值

（2）以下是一些常见的EPSG代码及其对应的坐标系：

1）EPSG:4326 - 这是WGS84地理坐标系的代码，目前是最流行的地理坐标系统之一。它广泛应用于GPS设备和在线地图服务中，因为GPS定位就是基于这个坐标系统的34。
2）EPSG:3857 - 网络墨卡托投影，适用于在线地图服务，如谷歌地图、百度地图和高德地图等。这种投影方式因其在互联网上的高效性而被广泛采用；
EPSG:4490 - 国家2000经纬度坐标系（CGCS2000）。

（3）epsg.io网站，以便用户可以轻松地查找特定EPSG代码对应的投影参数和其他相关信息:epsg代码查询网站

2-代码实现

2.1 数据介绍

2.2 代码实现

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

ds = gdal.Open(r"GDAL_testing_data\JAX_IMG1_PAN.TIF")

# ============================1.坐标系变换=================================
dsReprj = gdal.Warp("demReprj.tif", ds, dstSRS = "EPSG:4490")

# ============================2.影像重采样=================================
dsRes = gdal.Warp("demRes.tif", ds, xRes = 1, yRes = 1, 
                  resampleAlg = "bilinear")

# ============================3.影像裁剪=================================
# 确保栅格数据和 shapefile 具有相同的投影
dsClip = gdal.Warp("demClip.tif", ds, cutlineDSName = "star.shp",
                    cropToCutline = True, dstNodata = np.nan)

# ============================4.影像裁剪可视化=================================
array = dsClip.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(array)
plt.colorbar()
plt.show()

# close your datasets!
ds = dsClip = dsRes = dsReprj = None