python学习笔记1：python对遥感数据进行重投影

最新推荐文章于 2024-04-19 03:48:37 发布

Kkkary921

最新推荐文章于 2024-04-19 03:48:37 发布

阅读量411

点赞数 9

文章标签： python 学习笔记

本文链接：https://blog.csdn.net/Hey_everyday/article/details/134998052

版权

读取遥感影像数据并利用从参考影像的投影信息进行重投影，他们的空间范围一直

导入相关库

from osgeo import gdal
import numpy as np
import os
import fnmatch
import pandas as pd
import geopandas as gpd
import rasterio
from rasterstats import point_query

读写影像

def readData(envifile):
    envidata = gdal.Open(envifile)
    if envidata == None:
        print(envifile+' Open Failed')
    else:
        print(envifile+' Open Success')
    return envidata
def writeTiff(im_data, im_geotrans, im_proj, path):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float64

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape

    elif len(im_data.shape) == 2:
        im_data = np.array([im_data])
        im_bands, im_height, im_width = im_data.shape
    # 创建文件
    # 创建文件
    driver = gdal.GetDriverByName("Gtiff")
    dataset = driver.Create(path, int(im_width), int(
        im_height), int(im_bands), datatype)
    dataset.SetGeoTransform(im_geotrans)  # 写入放射变换函数
    dataset.SetProjection(im_proj)  # 写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])
    print(path + 'Write Success!')
    del dataset

重投影，根据一幅已有正确投影信息的影像对其他影像进行重投影

def reproject(inpath, outpath, referenceImg):

    filelist = []  # 用来存所有文件
    for root, dirs, files in os.walk(inpath):
        for name in files:
            if fnmatch.fnmatch(name, '*.tif'):
                filelist.append(root + '\\' + name)
    
    for imgfile in filelist:
        # imgfile = 'H:\\CHL\\new\\01.tif'
        dataset = readData(referenceImg)
        Img_width = dataset.RasterXSize  # 栅格矩阵的列数
        Img_height = dataset.RasterYSize
        Img_geotrans = dataset.GetGeoTransform()  # 获取放射矩阵信息
        Img_proj = dataset.GetProjection()  # 获取投影信息
        SavePath = '输出文件名.tif'
 
        godataset = readData(imgfile)
       data = godataset.ReadAsArray()
        data[data<-99] = 0  #把小于-99的值设置为0
        writeTiff(data, Img_geotrans, Img_proj, SavePath)

主函数

if __name__ == '__main__':
    inpath = r'M:\01'
    outpath = r'M:\11'
    referenceImg = r'reference.tif'#有正确投影信息的影像
    reproject(inpath, outpath, referenceImg)