如何将两个目录中的数组相乘以获得cloudfree多频带ras

2024-04-26 21:51:57 发布

您现在位置:Python中文网/ 问答频道 /正文

我想将目录中的4band光栅与相应的1波段云遮罩相乘,最终得到cloudfree 4band光栅。在我的尝试中,我首先将第一个波段乘以云遮罩。问题是脚本有时只返回一个文件,并显示我不理解的错误:

ERROR 1: Deleting D:/SR\20170207_4band.tif failed:
Permission denied



import pandas as pd
import gdal
import os
import numpy as np
import glob

SR_images = 'D:/SR'
sr = glob.glob(os.path.join(SR_images, '*.tif'))
reclass = 'D:/reclass'
mask = glob.glob(os.path.join(reclass, '*.tif'))
savedest = r'D:\maskedsr'


for e, i in zip(sr, mask):
        print(e)
        print(i)
        srraster = gdal.Open(os.path.join(SR_images, e))
        maskraster = gdal.Open(os.path.join(reclass, i))
        band1 = srraster.GetRasterBand(1).ReadAsArray()
        band2 = srraster.GetRasterBand(2).ReadAsArray()
        band3 = srraster.GetRasterBand(3).ReadAsArray()
        band4 = srraster.GetRasterBand(4).ReadAsArray()
        bandmask = maskraster.GetRasterBand(1).ReadAsArray()
        geotransform = srraster.GetGeoTransform()
        geoproj = srraster.GetProjection()
        xsize = srraster.RasterXSize
        ysize = srraster.RasterYSize
        format1 = "GTiff"
        driver = gdal.GetDriverByName(format1)
        dst_datatype = gdal.GDT_UInt16

        cal_band1 = band1 * bandmask
        name = e
        output_path = os.path.join(savedest, name)
        dst_ds = driver.Create(output_path, xsize, ysize, 1, dst_datatype)
        dst_ds.GetRasterBand(1).WriteArray(cal_band1)
        dst_ds.SetGeoTransform(geotransform)
        dst_ds.SetProjection(geoproj)
        dst_ds.FlushCache()
        dst_ds = None

谢谢


Tags: pathimportosdsglobdstjoingdal