py#gdal读写栅格

　1）栅格数读取

from osgeo import gdal
import numpy as np

#打开栅格数据集，只读
dataset=gdal.Open("E:/DATA/Image.TIF",gdal.GA_ReadOnly)
#读取图像信息
#数据格式
print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
                            dataset.GetDriver().LongName))
#栅格尺寸及波段数
print("Size is {} x {} x {}".format(dataset.RasterXSize,
                                    dataset.RasterYSize,
                                    dataset.RasterCount))
#投影坐标
print("Projection is {}".format(dataset.GetProjection()))
geotransform = dataset.GetGeoTransform()
if geotransform:
    print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
    print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))

#将栅格值转为array
#gdal中波段从1开始，不是从0开始
#读取之后array.shape=(RasterYSize,RasterXSize)，RasterY是纵轴方向即行数，RasterX为横轴方向即列数
band = dataset.GetRasterBand(1).ReadAsArray()

　2）栅格数据掩模处理

#提取小于等于1的影像，掩膜设为band>1
band_index=np.where(band>1)
band[band_index]=np.nan

　3）栅格数据写入

#先建立一个数据框架，确定输出数据的格式(GTiff)，文件名称，尺寸大小以及数值类型，尺寸大小先输入RasterX，即array的列数shape[1]，再输入RasterY，即array的行数shape[0]
gtiff_driver=gdal.GetDriverByName('GTiff')
out_ds=gtiff_driver.Create('test.tif',band.shape[1],band.shape[0],1,gdal.GDT_Float32)
#设置输出数据的坐标系
out_ds.SetProjection(dataset.GetProjection())
out_ds.SetGeoTransform(dataset.GetGeoTransform())
#将运算之后的波段信息写入数据框架中
out_ds.GetRasterBand(1).WriteArray(band)
#完成写入，清除数据集
out_ds=None

python小技巧

py#gdal读写栅格

相关

标签