1

I have images over 6 years of NDVI, NDWI, NDMI, MNDWI, SSI, IB and IB2 done with Sentinel 2. I'd like to have the mean value of each pixel over those years for each month and export it to be able to open it in QGis. For now my code looks like this :

# find the shape of an image

path= '/mnt/to/image.tif'
ds = gdal.Open(path)
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
pj = ds.GetProjection()

indices = ['NDVI', 'NDWI', 'NDWI', 'MNDWI', 'SSI', 'IB', 'IB2']
NAME_indices =['NDWI', 'NDWI', 'NDWI', 'MNDWI', 'SSI', 'IB', 'IB2']
NAME_mois =['janvier', 'fevrier', 'mars', 'avril', 'mai', 'juin', 'juillet', 'aout', 'sep', 'oct', 'nov', 'dec']

for y in indices:
    print(y, 'indice')
    for i in range(11):
        # df per month
        # get level values (0) because my df has several indexes
        print(NAME_mois[i])
        df_mois = df20[df20.index.get_level_values(0).month == i+1]
        # df per indice (and index lol)
        df_indices = df_mois.iloc[df_mois.index.get_level_values('Indice') == str(y)]
        
        # doing a column per pixel : this step took a long time 
        df_indices['Valeurs_clean'] = df_indices.loc[:,'Valeurs'].apply(lambda x: x.flatten()) 
        for z in range(len(df_indices.loc[:,"Valeurs_clean"][0])):
            df_indices["pixel" + str(z)] = df_indices['Valeurs_clean'].apply(lambda x:x[z])

I have 2 000 000 pixels per image and it took 400 minutes to have a column for 100 000 of them, then I create my map.

        # REFORMER LA CARTE 
        carte = np.reshape(df_indices, (height,width))
        print(carte)
        ### enregistrement image pour ensuite la réouvrir sur qgis 
        outFileName= "/mnt/Data/30_Stages_Encours/2023/ZonesHumides_Sarah/resultats/sentinel2/cartes/qgis/" + y + str(NAME_mois[i])+ ".tif"
        [rows, cols] = carte.shape
        driver = gdal.GetDriverByName("GTiff")
        outdata = driver.Create(outFileName, cols, rows, 1, gdal.GDT_Float32)
        print("Carte en cours de création")
        outdata.SetGeoTransform(gt) ##sets same geotransform as input
        outdata.SetProjection(pj) ##sets same projection as input
        outdata.WriteArray(carte)
        outdata.FlushCache() ##saves to disk!!
        outdata = None

I'd like to find a more efficient way to do this : mean of a pixel of the same column and list based on the index of the list maybe ? Do you have any idea. I'm also trying to put Dask into it but I'm new at this to it takes time.

To improve a little, I'm going to try to start with step before the rest to have to do it one time.

for z in range(len(df_indices.loc[:,"Valeurs_clean"][0])):
            df_indices["pixel" + str(z)] = df_indices['Valeurs_clean'].apply(lambda x:x[z])

I can't upload any image because it's too big.

sam
  • 47
  • 4

0 Answers0