GDAL修复Landsat ETM+影像条带
Landsat7 ETM+卫星影像由于卫星传感器故障,导致此后获取的影像出现了条带。如下图所示, 影像中均匀的布满条带。
使用GDAL修复影像条带的代码如下:
def gdal_repair(tif_name, out_name, bands): """ tif_name(string): 源影像名 out_name(string): 输出影像名 bands(integer): 影像波段数 """ # 打开影像文件 tif = gdal.Open(tif_name) # 根据文件类型获取对应的驱动程序 driver = gdal.GetDriverByName('GTiff') # 根据指定文件的驱动程序,使用现有数据集创建新的可写数据集 # 所有支持创建新文件的驱动程序都支持该`CreateCopy()`方法, # 但仅`Create()`部分支持该方法 # CreateCopy的第一个参数为目标文件名,第二个参数为源数据集 # 第三个参数的值是`0`或`1`,值是`0`。即使无法将原始数据准确地转换为目标数据,程序仍将执行 new_img = driver.CreateCopy(out_name, tif, 0) for i in tqdm(range(1, bands)): # 分别对每个波段处理 band = new_img.GetRasterBand(i) # 使用FillNodata对条带部分进行插值 gdal.FillNodata(targetBand = band, maskBand = band, maxSearchDist = 15, smoothingIterations=0) # 将修复好的波段写入新数据集中 new_img.GetRasterBand(i).WriteArray(band.ReadAsArray())
修复之后的效果图如下所示:
Opencv修复Landsat ETM+影像条带
使用opencv修复影像的代码如下:
def cv2_repair(tif_name): # 读取tif影像 tif_data = gdal_array.LoadFile(tif_name).astype('float32') # 获取掩膜 mask = tif_data.sum(axis=0) mask = (mask == 0).astype(np.uint8) bands = tif_data.shape[0] res = [] for i in tqdm(range(bands)): # cv.Inpaint(src, inpaintMask, dst, inpaintRadius, flags) # src:源图像,可以是8位、16位无符号整型和32位浮点型1通道或者8位无符号3通道 # inpaintMask:掩膜,8位无符号整型 # dst:和源图像具有一样大小的输出 # inpaintRadius:算法考虑的每个已修复点的圆形邻域的半径 # flags:修复算法类型,可选cv2.INPAINT_NS和cv2.INPAINT_TELEA repaired = cv2.inpaint(tif_data[i], mask, 3, flags=cv2.INPAINT_TELEA) res.append(repaired) return np.array(res)
修复之后的结果图:
使用opencv修复影像,速度要比Gdal慢许多,但修复质量更好。
Reference
https://www.bogotobogo.com/python/OpenCV_Python/python_opencv3_Image_reconstruction_Inpainting_Interpolation.php
https://gis.stackexchange.com/questions/151020/how-to-use-gdal-fillnodata-in-python