Python-gdal学习（1）栅格数据的操作

一、gdal简介

一个处理空间栅格数据的库

二、方法

0.引入gdal

from osgeo import gdal

1.注册驱动，即初始化对象

gdal.AllRegister() # 注册全部
driver=gdal.GetDriverByName("GTiff")  # 注册GTiff类型

2关于驱动

2.1查看支持的驱动数量

gdal.GetDriverCount()

2.2查看idx索引的驱动

driver=gdal.GetDriver(idx)
driver.ShortName
driver.LongName

3.读取遥感影像信息

3.1打开影像

src=gdal.Open('01.tif')

3.2查看对象可用操作

dir(src)

3.3获得栅格描述信息

src.GetDescription()

3.4获得栅格波段数

src.RasterCount #不带括号！

3.5获得数据宽高

width=src.RasterXSize
height=src.RasterYSize

3.6获取栅格六参数

src.GetGeoTransform()
#左上角x，水平分辨率，旋转，左上角y，旋转，竖直分辨率

3.7获取栅格投影

src.GetProjection()

3.8获取栅格元数据

src.GetMetadata()

3.9获取波段，从1开始索引

band1=src.GetRasterBand(1)

3.10查看波段基本信息：属性和方法

dir(band1)

3.10.1查看波段数据类型

band1.DataType

3.11波段宽高和空间范围.

band1.XSize,band1.YSize,band1.DataType

3.12波段数据的属性

band1.GetNoDataValue()
band1.GetMaximum()
band1.GetMinimum()
band1.ComputeRasterMinMax()

None
None
None
(0.0,255.0)

4.获取数据

4.1引入及查看常量

from osgeo import gdalconst
dir(gdalconst)

4.2访问数据集的数据，6个参数：左上角x,y，在x,y上分别取得像元数，重采样像元数，默认同参数3，4；索引时不要越界

src.ReadAsArray(2500,2500,3,3)
src.ReadRaster(2500,2500,3,3)

5.复制tif

new='01_copy.tif'

5.1方法一得到01_copy.tif

new_ds=driver.CreateCopy(new,src,0)

5.1.1加参数！

# 压缩方式，上面文件我创建失败
new_ds=driver.CreateCopy(dist,src,0,["TITLED=YES","COMPRESS=PACKBITS"])

new_ds=driver.CreateCopy(dist,src,0,["TITLED=NO","COMPRESS=PACKBITS"])

5.2方法二按像元存储顺序

# 像元存储
new_ds=driver.CreateCopy(new,src,0,["INTERLEAVE=BAND"])
# 按不同方式存储：BSQ、BIP
new_ds=driver.CreateCopy(new,src,0,["INTERLEAVE=PIXEL"])

6创建tif

6.1创建单波段tif

new_ds=driver.Create(new,512,512,1,gdal.GDT_Byte) #创建单波段空tif

6.1.2空间参数信息单独处理。添加空间范围和坐标系

import numpy,osr
new_ds.SetGeoTransform([444720,30,0,3751320,0,-30]) #安一个六参数
srs=osr.SpatialReference()
srs.SetUTM(11,1) # 设置投影带
srs.SetWellKnownGeoCS('NAD27') #空间参考
dst_ds.setProjection(srs.ExportToWkt())
raster=numpy.zeros((512,512))
new_ds.GetRasterBand(1).WriteArray(raster)

6.2.1直接创建多波段tif

new_ds=driver.Create(new,512,512,3,options=["INTERLEAVE=PIXEL"])
#写入数据
new_ds.WriterRaster(0,0,512,512,datas.tostring(),512,512,band_list=[1,2,3])

new_ds.FlushCache() #刷新缓存

6.2.2分波段处理

datas=[]
for i in range(3):
    band=dataset.GetRasterBand(i+1)
    data=band.ReadAsArray(0,0,width,height)
    datas.append(numpy.reshape(data,(1,-1))) #拼接
   datas=numpy.concatenate(datas)

6.3输出影像

new=driver.Create('/wenjain.tif',width,height,3,option=["INTERLEAVE=PIXEL"])
new.WriteRaster(0,0,width,height,datas.tostring(),width,height,band_list=[1,2,3])
# 用ReadAsArray读数据，需要用tostring转回数据。若使用ReadRaster则不需要
new.FlushCache()

6.4影像的空间参数信息单独处理。例如导入一个NAD27的空间参考

new_ds.SetGeoTransform(src.GetGeoTransform)
srs=osr.SpatialReference()
srs.SetUTM(11,1)
srs.SetWellKnownGeogCE('NAD27')
new_SetProjection(srs.ExportToWkt())

6.5建立影像金字塔

gdal.SetConfigOption('HFA_USE_RRD','YES')
# 建立金字塔
new.BuildOverviews(overviewlist=[2,4,8,16,32,64,128])