一、读取

安装方法:

  1. 手动下载whl安装包https://www.lfd.uci.edu/~gohlke/pythonlibs/
  2. 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文件
Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐