python gdal读写tif栅格数据
参考:https://www.cnblogs.com/ninicwang/p/11533066.html# -*- coding: utf-8 -*-from osgeo import gdaldataset = gdal.Open(r'F:\MOD11A1.A2018152.LST_Day_1km.tif')#栅格矩阵的列数im_width = dataset.RasterXSize#栅格矩阵的
·
一、读取
安装方法:
- 手动下载whl安装包https://www.lfd.uci.edu/~gohlke/pythonlibs/
- pip install 安装包
# -*- coding: utf-8 -*-
from osgeo import gdal
dataset = gdal.Open(r'F:\MOD11A1.A2018152.LST_Day_1km.tif')
#栅格矩阵的列数
im_width = dataset.RasterXSize
#栅格矩阵的行数
im_height = dataset.RasterYSize
#波段数
im_bands = dataset.RasterCount
#仿射矩阵,左上角像素的大地坐标和像素分辨率。
#共有六个参数,分表代表左上角x坐标;东西方向上图像的分辨率;如果北边朝上,地图的旋转角度,0表示图像的行与x轴平行;左上角y坐标;
#如果北边朝上,地图的旋转角度,0表示图像的列与y轴平行;南北方向上地图的分辨率。
im_geotrans = dataset.GetGeoTransform()
#地图投影信息
im_proj = dataset.GetProjection()
#读取某一像素点的值
#(1)读取一个波段,其参数为波段的索引号,波段索引号从1开始(我打开的这幅图像只有一个波段)
band = dataset.GetRasterBand(1)
#(2)用ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>),读出从(xoff,yoff)开始,大小为(xsize,ysize)的矩阵。以下为读取整幅图像
im_datas = band.ReadAsArray(0,0,im_width,im_height)
#释放内存。如果不释放,在arcgis或envi中打开该图像时显示文件已被占用
del dataset
二、写入
方法一:参照已有tif数据写
from osgeo import gdal
# 打开原始TIFF文件
input_ds = gdal.Open('original.tif', gdal.GA_ReadOnly)
# 获取原始文件的GeoTransform和Projection
geotransform = input_ds.GetGeoTransform()
projection = input_ds.GetProjection()
# 假设我们要创建一个新的TIFF文件,它具有与原始文件相同的尺寸和波段数
# 注意:这里我们只是创建了一个空的数据集,你可能需要填充数据
driver = gdal.GetDriverByName('GTiff')
width = input_ds.RasterXSize
height = input_ds.RasterYSize
bands = input_ds.RasterCount
data_type = input_ds.GetRasterBand(1).DataType
# 创建新的数据集
output_ds = driver.Create('new.tif', width, height, bands, data_type)
# 设置新的数据集具有与原始数据集相同的GeoTransform和Projection
output_ds.SetGeoTransform(geotransform)
output_ds.SetProjection(projection)
# ... 在这里你可以填充数据到output_ds ...
output_ds.GetRasterBand(1).WriteArray(data)
# 关闭数据集
del input_ds # 或者 input_ds = None
del output_ds # 或者 output_ds = None
方法二:无参照tif文件写
from osgeo import gdal,osr
def output_tif(output_path, lon, lat, lon_res, lat_res, arr_2d):#lon,lat均为一维数据,且无重复数据
driver = gdal.GetDriverByName('GTiff')
out_tif = driver.Create(output_path, len(lon), len(lat), 1, gdal.GDT_Float32) # 设置影像的显示范围#-Lat_Res一定要是-的
geotransform = (np.min(lon), lon_res, 0, np.max(lat), 0, -lat_res)
out_tif.SetGeoTransform(geotransform) # 获取地理坐标系统信息,用于选取需要的地理坐标系统
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84",AUTHORITY["EPSG","4326"]
out_tif.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息#数据写出
out_tif.GetRasterBand(1).WriteArray(arr_2d) # 将数据写入内存,此时没有写入硬盘
out_tif.FlushCache() # 将数据写入硬盘
out_tif = None # 注意必须关闭tif文件
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
已为社区贡献5条内容
所有评论(0)