Gdal读写栅格数据

Gdal读写栅格数据

仅提供读取数据和写入数据功能，后续可能更新

from osgeo import gdal
# py2 import gdal


class Gdal:
    def __init__(self, filepath):
        self.filepath = filepath
        self.im_width = None
        self.im_height = None
        self.im_data = None
        self.im_geotrans = None
        # 假设指北针向上的影像，GT2和GT4参数是0，而GT1是象元宽，GT5是象元高，
        # （GT0，GT3）点位置是影像的左上角。
        self.im_proj = None  # 存储投影信息

    def readtif(self):
        dataset = gdal.Open(self.filepath)
        self.im_width = dataset.RasterXSize
        self.im_height = dataset.RasterYSize
        self.im_proj = dataset.GetProjection()
        self.im_geotrans = dataset.GetGeoTransform()
        self.im_data = dataset.ReadAsArray(0, 0, self.im_width, self.im_height)
        del dataset
        return self.im_width, self.im_height, self.im_data.astype(float)

    def writeTif(self, newpath, datatype, im_data, select, nodata=None):  # datatpye 如无符号16为整型
        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)
        driver = gdal.GetDriverByName('GTiff')
        # 多波段计算结果以单波段形式保存
        if select == 1:
            im_bands = 1
        new_dataset = driver.Create(newpath, im_width, im_height, im_bands, datatype)
        new_dataset.SetGeoTransform(self.im_geotrans)
        new_dataset.SetProjection(self.im_proj)
        band = new_dataset.GetRasterBand(1)
        band.SetNoDataValue(nodata)
        band.FlushCache()
        if im_bands == 1:
            new_dataset.GetRasterBand(1).WriteArray(im_data)
        else:
            for band in range(1,im_bands+1):
                new_dataset.GetRasterBand(band).WriteArray(im_data)
        del new_dataset