程序目的:python语言,使用gdal读取哨兵2影像,转换并写入为geotiff格式。由于哨兵存储格式复杂,使用循环语句按子数据集读取波段并写入新影像中。哨兵的数据结构为[2, 3, 4, 8], [5, 6, 7, 10, 11, 12], [1, 9]三个子集,各自集分辨率和影像大小不同。
错误:新写入的tif影像数据不全,且写入的波段也为无法解译的噪点。通过查看,第一子集的数据成功保存,但为噪点无法解译;第二三子集为空。
报错:pycharm编译时报三类错误:
ERROR 1:Cannot find proj.db #虽然报错,投影似乎还是正确的
ERROR 5: Destination band 1 appears to be read-only.
ERROR 1: Buffer too small
代码:
# 将波段进行提取,并存储至Geotiff文件
def ReadWrite(fileset,lenrt):
# 新建一个2维数组,用以给写入的波段排序(使用numpy)
bandlist = np.array([[2, 3, 4, 8], [5, 6, 7, 10, 11, 12], [1, 9]], dtype=np.object)# 10实际为8a
# numpy写入的数组不能直接定位,因为其组成为3个list
# 循环读取并写入数组
for i in range(lenrt): # 1循环打开影像并读取所有子数据集,除全色以外
fullset = gdal.Open(fileset[i]) # 注意是i
subset = fullset.GetSubDatasets()
IN = 'E:/process/Ubuntu/' + str(i+1) + '.tif'
cankao = 'E:/process/Ubuntu/B1.img'
ImageE = CreateE(cankao, IN)
print('打开第'+str(i+1)+'幅影像,路径为'+str(fileset[i]))
for j in range(0, len(subset) - 1): # 2循环打开单一影像所有子数据集
print('打开第'+str(j+1)+'个子数据集,数据集包括波段'+str(bandlist[j]))
R10m = gdal.Open(subset[j][0])
width = R10m.RasterXSize # 读取X
height = R10m.RasterYSize # 读取Y
if width==10980 and height==10980:
b10m = R10m.ReadAsArray(0, 0, width, height) # 将子数据集读取为ndarray&重采样
else:
R10n = R10m
gdal.Warp(R10n, R10m, width=10980,height=10980,
resampleAlg=gdal.GRIORA_Bilinear)
width = R10n.RasterXSize # 读取X
height = R10n.RasterYSize # 读取Y
b10m = R10n.ReadAsArray(0, 0, width, height)
np.reshape(b10m, (-1, 1)) # 归一化
b10m = b10m.tobytes()
subsetlist = bandlist[j] # 因为bandlist是numpyarray格式,所以要预先取出其中的list
# 保存所有波段为一个TIFF
ImageE.WriteRaster(0, 0, width, height, b10m, buf_xsize=10980, buf_ysize=10980, band_list=subsetlist)
# WriteRaster是与Readraster相对应的函数,它将二进制形式的波段数据写入空数据集中。WriteRaster(
# self, xoff, yoff, xsize, ysize,
# buf_string, buf_xsize, buf_ysize, buf_type,
# band_list, buf_pixel_space, buf_line_space, buf_band_space)
# x y off为数据的偏移;x y size是创建波段的大小;
j = j + 1
if j == len(subset)-1: # 2读取结束
break
ImageE.FlushCache()
print('影像'+str(fileset[i])+'转换成功')
i = i + 1
if i == lenrt - 2: # 1读取结束
break
return 0
如果查看不便,请私信我发送代码文件gei'ni