python批量转换:未知类型二进制净生产力(NEP)打开遥感定量反演产品数据,转换为tif,可用于通用GIS、打开遥感软件。
完整实例下载:批量分析未知遥感产品数据并转换为tif
问题描述:遥感没有头文件,不能使用envi软件直接读取时,不能使用gdal读取。图像格式未知,存储模式未知(只知行列、字节数、二进制无符号短整形)。
数据集从 1981 年 1 月 1 日开始至 2019 年 12 月 31 日,每年 365期数据。为了节省存储空间,提高下载速度,压缩数据一年 365 期数据压缩成文件。如2019NEP.zip 包含在文件中 365 个文件。文件采用 NEP_YYYY_ddd.img 格式命名,其中 YYYY 表示年 份,ddd 表示日序。如 NEP_2019_1.img 是 2019 年第 1 天的 NEP;NEP_2019_101.img 是 2019 年累积到 101 天的 NEP,比例因子为 0.1。减少相邻两天的数据,可以得到每天NEP值。如将从NEP_2019_101.img从文件读出的值减去 NEP_2019_100.img 文件读取的值乘以比例因子是第一个 101 天的 NEP(g C m-2 d-1)。NEP 正值表示生态系统吸收大气 CO2,NEP 负值表示生态系统向大气排放 CO2。 数据存储模式:有符号短整形二进制文件,每个大小为 2 字节/格点×4950 格点/行×2090 行=20691000 字节。 数据存储顺序如下:BSQ
解决方案的整个过程: (1)确定每个像素点所占字节数: 读取以下代码:图像中的总字节数,运行结果为sum。由于一幅影像是由sumPixel=“行 * 列 * 波段数” 因此,推断每个像素点都是由像素点组成的"sum/sumPixel"字节组成。
# -*- coding: utf-8 -*- import numpy as np def test_readbil(): dir =r"C:\Users\HP\Desktop\imgFolder\NEP_2019_364.img" f = open(dir,'rb') n = 0 while True: s = f.read(1) if s == b"": break n = n 1 print(n) #输出:42105600
(2)读取并转存为tif格式图像
根据以下三种存储格式进行测试: 一、BIL(Band Interleaved by Line):根据存储一行的所有波段,然后存储下一行的存储方法。 二、BSQ格式(Band Sequential Format):存储下一波段整幅图像后,按存储单波段整幅图像的方法存储。(一般为此方法)。 三、BIP格式(Band Interleaved by Pixel):存储一个像素的所有波段,然后存储下一个像素的所有波段。
import numpy as np import gdal import os # 依据BSQ存储规则,根据存储单波段整个图像后存储下一波段的存储方法提取并存储在数组中。 def read_as_bsq(dataarr,bands, rows, col): imgarr = np.zeros((bands,rows,col)) for b in range(bands): #取出每个波段 start = b * rows * col end = start rows * col arrtem = dataarr[start:end]
for r in range(rows): #一维数组按行取出,存入对应三维数组。
start2 = r * col
end2 = start2 + col
imgarr[b,r,:] = arrtem[start2:end2]
return imgarr
def test_writetotif(img_dir,tif_savedir):
#读取二进制数据并转换成Int16类型的数组
#img_dir = r"C:\Users\HP\Desktop\imgFolder\NEP_2019_364.img"
f = open(img_dir, 'rb')
#数据提取关键代码:
省略(联系我)
#将提取的数组存储为tif格式图像.
#注意这里未设置地理坐标和仿射变换信息,所以不能用ENVI等软件打开。
#tif_savedir = r"C:\Users\HP\Desktop\tifFolder\NEP_2019_364.tif"
datatype = gdal.GDT_Int16 #gdal.GDT_UInt16
bands, high, width = imgarr.shape
driver = gdal.GetDriverByName("GTiff")
datas = driver.Create(tif_savedir, col, rows, bands, datatype)
#设置地理坐标和仿射变换信息
# datas.SetGeoTransform(image_geotrans)
# datas.SetProjection(image_projetion)
for i in range(bands):
datas.GetRasterBand(i + 1).WriteArray(imgarr[i])
del datas
#print("save succfully")
if __name__ =="__main__":
test_writetotif(imgPath, tifPath)
如果有疑问或者需要相关批量的遥感功能,请私信联系我。