赞
踩
想了一下,先说一下读写tif这两个基本操作吧
1、读取tif为矩阵
#获取文件句柄
handle = gdal.Open(‘123.tif’)
#获取文件的列数、行数和波段数
col = handle.RasterXsize
row = handle.RasterYsize
#获取放射信息和投影信息
geotrans = handle.GetGeoTransform()
geoinfo = handle.GetProjectionRef()
#获取数据,其中xoff,yoff是想要读取的部分原点相对于源数据远点的位置(xoff是列偏移量,yoff是行偏移量,一般设置成0,0)col,rowsize是要读取数据的列数和行数
dataset = handle.ReadAsArray(xoff,yoff,xsize,rowsize) #一般使用dataset = handle.ReadAsArray(0,0,col,row)
#上面获得的是三维矩阵,也可以先获得某个波段(波段号从1开始),这里面
bandset = handle.GetRasterBand(bandindex)
dataset = bandset.ReadAsArray(xoff,yoff,xsize,rowsize)
2、矩阵写入tif
#注册一个driver
driver = gdal.GetDriverByName('GTiff')
#创建数据集,参数为文件名,列数,行数,波段数,数据类型
dataset = driver.Create(filename,width,height,band,datatype)
#设置放射矩阵,leftlon,leftlat 分别是左上角的经度纬度,lonres,latres分别是经度纬度的分辨率,剩下两个0的参数是图像的旋转角度,对于北方朝上的图像来说,这两个值一般是0
dataset.SetGeoTransform(leftlon,lonres,0,leftlat,0,latres)
#设置投影信息,参数是OpenGis的 WKT字符串格式
dataset.SetProjection(WKT)
#写数据,bandindex是波段序号,一般从1开始,data是矩阵数据
dataset.GetRasterBand(bandindex).WriteArray(data)
3、写成tif时出现的数值错误的问题
之前我部门头儿写了个python使用gdal库读写的tif的类,给我使用,发现了他把所有数据转换成了int型,然后写入,当时还想着说估计是疏忽了,所以给他加了一条判断,代码如下
跑完程序结果如下
分隔栏上面的是源数据,下面的事写入dataset之后的样子,没错,两边都是float怎么会这样呢,按理说是不存在类型转换之间的问题。我还查看了一下
分隔栏上面是源数据类型,底下是写入数据的类型,我查看了一下f4这个类型的name就是float32。没道理,所以就尝试着将转换的矩阵转换成了numpy类型,像这样。
得到的结果,就对了,看来平时的float和numpy的float是不一样的,至少,在gdal的内部转化中是不一样的,也可能是gdal的这个ReadAsArray对python自带的数据类型并没有兼容。
之前发现了这个问题就放一边了,没时间想咋回事。要不是碰到今天这个问题数据,不知道还要搁置多久。。。
最后求关注,求点赞,欢迎大家关注我的公众号
记录所学所用,包括但不限于遥感、地信、气象、生态环境,机器学习知识,相关文献阅读,编程代码实现。偶尔荒腔走板的聊聊其他。欢迎不同领域的朋友们加入进来,多多交流。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。