I want to automate this code below for other 99 images, the problem that I have is that each image is in its specific folder as it's mentioned in the picture below.

enter image description here

the images are very big and I want to clip them in a specific area.

Is that possible to make it automatic also ?

    import rasterio
import os
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np

os.chdir(r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP')

band4 = rasterio.open(r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP\clip_RT_L1C_T29SPR_A003995_20160328T111353_B04.tif') #red
band5 = rasterio.open(r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP\clip_RT_L1C_T29SPR_A003995_20160328T111353_B08.tif') #nir


red = band4.read(1).astype('float32')
nir = band5.read(1).astype('float32')


plt.colorbar (ndvi) #I have also a problem with colorbar as it didn't show properly

Something like this:

path = r'Z:\New folder\S2\L1C_T29SPR_A003995_20160328T111353_2016-03-28\New folder\CLIP\'
dirContents = os.listdir(path)

for file in dirContents:
    if os.path.isdir(file):
        subDir = os.listdir(file)
        # Assuming only two files in each subdirectory, bands 4 and 8 here
        if "B04" in subDir[0]:
            band4 = rasterio.open(subDir[0])
            band5 = rasterio.open(subDir[1])
            band4 = rasterio.open(subDir[1])
            band5 = rasterio.open(subDir[0])

        red = band4.read(1).astype('float32')
        nir = band5.read(1).astype('float32')

        # And do the rest of your analysis here. NDVI calculation, clipping, etc.
  • Thank you for your reply. Can you be more precise, I'm very new in python. – M.Yassine Mar 27 '19 at 18:38

I would recommend using rioxarray for this task.

Here is how I would do it:

from pathlib import Path
import rioxarray

base_folder = Path("Z:\\New folder\\S2\\")

all_folders = base_folder.glob("*")

geometries = [
        'type': 'Polygon',
        'coordinates': [[
            [425499.18381405267, 4615331.540546387],
            [425499.18381405267, 4615478.540546387],
            [425526.18381405267, 4615478.540546387],
            [425526.18381405267, 4615331.540546387],
            [425499.18381405267, 4615331.540546387]

def load_clipped(in_path, in_geometries):
    in_data = rioxarray.open_rasterio(in_path).sel(band=1).astype('float32')
    return in_data.rio.clip(in_geometries)

for folder in all_folders:
    if not folder.is_dir():
        red_path = folder.rglob("*_B04.tif")[0]
        base_dir = red_path.parent
        red = load_clipped(red_path, geometries)
        nir = load_clipped(base_dir.glob("*_B08.tif")[0], geometries)
    except IndexError:

    ndvi = ((nir-red)/(nir+red)).where((nir+red)!=0., 0)
    ndvi.attrs = red.attrs
    ndvi = ndvi.rio.write_crs(red.rio.crs)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.