||
(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)
【参考】
https://www.jb51.net/article/151644.htm
点滴分享,福泽你我!Add oil!
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 14:39
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社