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 |
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 【.NET】调用本地 Deepseek 模型
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库