NJU1healer的个人博客分享 http://blog.sciencenet.cn/u/NJU1healer

博文

利用python GDAL库读写geotiff格式的遥感影像方法

已有 6497 次阅读 2020-9-15 17:10 |个人分类:Python|系统分类:科研笔记

(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!



https://wap.sciencenet.cn/blog-3428464-1250704.html

上一篇:Python代码调试错误集合(1)
下一篇:查询Anaconda安装路径、安装包位置及授予访问权限
收藏 IP: 211.162.81.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

全部作者的其他最新博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 17:12

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部