(1)利用python GDAL库读写geotiff格式的遥感影像方法，具有很好的参考价值，不错！

from osgeo import gdal

import numpy as np

def read_tiff(inpath):

ds=gdal.Open(inpath)

row=ds.RasterXSize

col=ds.RasterYSize

band=ds.RasterCount

geoTransform=ds.GetTransform()

proj=ds.GetTransform()

data=np.zeros([row,col,band])

for i in range(band):

dt=ds.GetRasterBand(1)

data[:,:,i]=dt.ReadAsArray(0,0,col,row)

return data,geoTransform,proj

def array2raster(outpath,array,geoTransform,proj):

cols=array.shape[1]

rows=array.shape[0]

driver=gdal.GetDriverByName('Gtiff')

outRaster=driver.Create(newRasterfn,cols,rows,1,gdal.GDT_Byte)

outRaster.SetGeoTransform(geoTransform)#参数2,6为水平垂直分辨率，参数3,5表示图片是指北的

outband=outRaster.GetRasterBand(1)

outband.WriteArray(array)

outRaster.SetProjection(proj)#将几何对象的数据导出为wkt格式

outRaster.FlushCache()

if __name__=="_main_":

data,geoTransform,proj=read_tiff('d:/a.tif')

array2raster("d:/b.tif",np.zeros[2400,2400],geoTransform,proj)

(2)具体细节参数

#1*************波段组成图像，波段指挥颜色************

import gdal  #导入库

dataset=gdal.open("filename")  #文件名，如*.tif

dir(dataset)     #可通过dir()来显示信息，或者某个函数的用法or  print(dataset)

#基本函数

dataset.GetDescription()   #获得栅格的描述信息

dataset.RasterCount        #获得栅格数据集的波段数

band=dataset.GetRasterBand(1)  #获得栅格数据集的波段

width=dataset.RasterXSize  #读取图像的宽度，x方向上的像素个数

height=dataset.RasterYSize  #读取图像的高度，y方向上的像素个数

#读取二进制

dataset.ReadRaster(xoff,yoff,xsize,ysize,buf_xsize=None,buf_ysize=None,buf_type=None,band_list=None)

#xoff,yoff指定想要读取的部分原点位置在整张图像中距离全图原点的位置

#xsize和ysize指定要读取部分图像的矩形大小

#buf_xsize,buf_ysize代表缩放大小，buf_type可设置读取的数据类型

#读取数组

dateset.ReadAsArray(xoff,yoff,xsize,ysize)

#举个例dataset.ReadRaster(230,270,10,10)#把图像中位于230，270，宽度为10高度为10的数据读取出来了

band.XSize

band.YSize  #波段图像的宽和高

band.DataType  #图像中实际数值的数据类型，具体的数据类型定义在gdalconst模块里，需import gdalconst

band.GetNoDataValue() #获取无意义值

band.GetMaximum()

band.GetMininum()  #由于是文件格式没有固有的最值，所以不会显示出来，用下面这个函数

band.ComputeRasterMinMax()#能算出排除了无意义值后的最值

band.GetRasterColorInterpretation()

gdalconst.GCI_PaletteIndex

colormap=band.GetRasterColorTable() #获得颜色表

colormap.GetCount()   #获得颜色的数量

colormap.GetPaletteInterpretation() #知道获得的颜色表是什么颜色表，比如是RGB颜色表

#比如GPI_GRAY,GPI_RGB,GPI_CYMP,GPI_HLS

gdal.GPI_RGB

for i in range(colormap.GetCount()):

print(colormap.GetColorEntry(i))  #获得颜色表中的值，有四个值，其中前三个有意义

#2**************细节**************

band.ReadAsArray(xoff,yoff,win_xsize,win_ysize,buf_xsize,buf_ysize)

#xoff,yoff是取值窗口的左上角在实际数据中所处象元的xy位置。

#win_xsize,win_ysize是取值窗口覆盖的区域大小

#buf_xsize,buf_ysize是取值窗口取出数组进行缩放后数组的大小，缩的时候是取周围点的平均值，如果不设置，则跟3、4参数一致

#例子1：横向读取和纵向读取的效率

import gdal

import time

dataset=gdal.Open("filename")

band=dataset.GetRasterBand(1)

width=dataset.RasterXSize

height=dataset.RasterYSize

bw=128

bh=128

bxsize=width/128

bysize=height/128

start=time.time()

band.ReadAsArray(0,0,width,height)

print(time.time()-start)

start2=time.time()

for i in range(bysize):

for j in range(bxsize):

band.ReadAsArray(bw*j,bh*i,bw,bh)

print(time.time()-start2)

#调换循环顺序

for j in range(bxsize):

for i in range(bysize):

band.ReadAsArray(bw*j,bh*i,bw,bh)

【参考】

点滴分享，福泽你我！Add oil！

转载本文请联系原作者获取授权，同时请注明本文来自张伟科学网博客。

链接地址：http://blog.sciencenet.cn/blog-3428464-1250704.html

上一篇：Python代码调试错误集合(1)

下一篇：查询Anaconda安装路径、安装包位置及授予访问权限