py#gdal读写栅格

 1)栅格数读取

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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
2
3
#提取小于等于1的影像,掩膜设为band>1
band_index=np.where(band>1)
band[band_index]=np.nan

 3)栅格数据写入

1
2
3
4
5
6
7
8
9
10
#先建立一个数据框架,确定输出数据的格式(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

  

posted @   MisakaYier  阅读(184)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 【.NET】调用本地 Deepseek 模型
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
点击右上角即可分享
微信分享提示