[Python知识库]python使用gdal读写栅格数据

import os
import gdal
import numpy as np
import pandas as pd

def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName + "文件无法打开")
        
    width = dataset.RasterXSize 
    height = dataset.RasterYSize 
    bands = dataset.RasterCount 

    data = dataset.ReadAsArray()
    geotrans = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    return width, height, bands, data, geotrans, proj

def writeTiff(im_data, im_geotrans, im_proj, path):
    im_data[np.isinf(im_data)] = 0  #处理异常值
    im_data[np.isneginf(im_data)] = 0
    im_data[np.isnan(im_data)] = 0
    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_Float32
    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)
    if(dataset!= None):
        dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
        dataset.SetProjection(im_proj) #写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset

if __name__ == '__main__':
    fileName = r'F:\package\xnd\index\sy_index_result.tif'
    path = r'F:\package\xnd\index\sy_index_result_pro.tif'
    width, height, bands, data, geotrans, proj = readTif(fileName)
    writeTiff(data,geotrans, proj, path)