编程_栅格影像去除投影

ytkz2023-12-052024-09-09

地理坐标系一般是指由经度、纬度和相对高度组成的坐标系，能够标示地球上的任何一个位置。经度和纬度常合称为经纬度，把球面上的经纬度显示在平面地图上需要采用某种地图投影

地理坐标系

由于地面高低的不同和地球形状的不正规，天体测量所得的信息不足用于明确地计算地理位置。一般采用某个标准的大地测量系统所规定的地理坐标系统，目前最常用的标准是WGS84，它被卫星导航系统如美国GPS、中国的北斗导航系统使用。

投影坐标系

投影坐标系最出名的墨卡托投影，这里不过多讨论。

需求

现在有一份栅格影像tif文件，需要把这份栅格文件去除投影坐标，变为地理坐标。准确地说是，让这景tif的坐标系从GCS2000 变为WGS84 。

代码

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2023/12/5 14:40 
# @File : reprojection.py


# 栅格数据可以重投影，但比矢量数据投影更复杂。
# 栅格数据需要处理栅格数据中像元会弯曲和移动的事实，一对一的映射并不存在
# 通常用最近邻域插值法、双线性插值和三次卷积插值法进行插值
# 方法：1.gdalwarp 2.AutoCreateWarpedVRT

from osgeo import gdal, osr
import os

def reproject(infile, outfile):
    if os.path.exists(os.path.dirname(outfile)):
        pass
    else:
        os.makedirs(os.path.dirname(outfile))

    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS('WGS84')  # UTM转无投影，即地理坐标系
    old_ds = gdal.Open(infile)
    vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(),
                                      gdal.GRA_Bilinear)  # 第一个默认值None，使用源栅格数据的srs；第二个如果None，表示不发生重投影
    gdal.GetDriverByName('gtiff').CreateCopy(outfile, vrt_ds)  # 该函数返回一个数据集对象，用CreateCopy函数保存为gtiff

if __name__ == '__main__':

    file = r'D:\img1.tif'
    outfile = r'D:\img1_wgs84.tif'
    reproject(file, outfile)