遥感+python 1.5 重投影

原创已于 2023-08-02 15:03:32 修改 · 464 阅读

3 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#python

于 2023-07-11 17:31:01 首次发布

遥感+python 专栏收录该内容

7 篇文章

订阅专栏

本文介绍了遥感影像处理中的重投影概念，通过Python使用GDAL库来读取、处理和输出影像，并详细阐述了如何进行坐标转换，包括设置新投影、创建转换对象以及更新地理变换参数。

遥感+python 1.5 重投影

一、重投影概念

当考虑两幅遥感影像、矢量地图等的坐标信息时，我们需要考虑其所在的投影坐标系。若两投影坐标系不相同则需要进行重投影变换。

地理坐标系、大地坐标系与地图投影与重投影详解

二、代码实现

2.1 读取影像

利用gdal读取遥感影像

from osgeo import gdal

def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    if dataset == None:
        raise Exception(f"cant find/open {filename}")
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_Band = dataset.RasterCount  # 栅格矩阵的波段数

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵

    del dataset  # 关闭对象，文件dataset
    return im_proj, im_geotrans, im_data, im_width, im_height, im_Band

2.2 输出影像

DType2GDAL = {"uint8": gdal.GDT_Byte,
              "uint16": gdal.GDT_UInt16,
              "int16": gdal.GDT_Int16,
              "uint32": gdal.GDT_UInt32,
              "int32": gdal.GDT_Int32,
              "float32": gdal.GDT_Float32,
              "float64": gdal.GDT_Float64,
              "cint16": gdal.GDT_CInt16,
              "cint32": gdal.GDT_CInt32,
              "cfloat32": gdal.GDT_CFloat32,
              "cfloat64": gdal.GDT_CFloat64}


def write_img(filename, im_proj, im_geotrans, im_data):
    # 判断栅格数据的数据类型
    if im_data.dtype.name in DType2GDAL:
        datatype = DType2GDAL[im_data.dtype.name]
    else:
        datatype = gdal.GDT_Float32
    # 判读数组维数
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape
    # 创建文件
    if not pathlib.Path(filename).parent.exists():
        pathlib.Path(filename).parent.mkdir(parents=True, exist_ok=True)
    driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset

2.3 坐标转换

class GEOchange(object):

    def __init__(self, toEPSG):
        self.EPSG = toEPSG
        self.to_crs = osr.SpatialReference()
        self.to_crs.ImportFromEPSG(toEPSG)

    def run(self, infile, outfile):
        im_proj, im_geotrans, im_data, self.im_width, self.im_height, self.im_Band = read_img(infile)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(im_proj)
        self.Transformation = osr.CoordinateTransformation(srs, self.to_crs)
        geotrans = self.setGeotrans(im_geotrans)
        write_img(outfile, self.to_crs.ExportToWkt(), geotrans, im_data)

    def setGeotrans(self, im_geotrans):
        lon, lat = self.imagexy2geo(im_geotrans, 0, 0)
        coords00 = self.Transformation.TransformPoint(lat, lon)
        lon, lat = self.imagexy2geo(im_geotrans, self.im_height, 0)
        coords01 = self.Transformation.TransformPoint(lat, lon)
        lon, lat = self.imagexy2geo(im_geotrans, 0, self.im_width)
        coords10 = self.Transformation.TransformPoint(lat, lon)

        trans = [0 for i in range(6)]
        trans[0] = coords00[0]
        trans[3] = coords00[1]
        trans[2] = (coords01[0] - trans[0]) / self.im_height
        trans[5] = (coords01[1] - trans[3]) / self.im_height
        trans[1] = (coords10[0] - trans[0]) / self.im_width
        trans[4] = (coords10[1] - trans[3]) / self.im_width
        return trans

    def imagexy2geo(self, im_geotrans, row, col):
        px = im_geotrans[0] + col * im_geotrans[1] + row * im_geotrans[2]
        py = im_geotrans[3] + col * im_geotrans[4] + row * im_geotrans[5]
        return px, py

<返回目录