1

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

plot.show(band5)
band4.dtypes

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

ndvi=np.where(
    (nir+red)==0.,
    0,
    (nir-red)/(nir+red))
ndvi[:5,:5]

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

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])
        else:
            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 at 18:38
0

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():
        continue
    try:
        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:
        continue

    ndvi = ((nir-red)/(nir+red)).where((nir+red)!=0., 0)
    ndvi.attrs = red.attrs
    ndvi = ndvi.rio.write_crs(red.rio.crs)
    ndvi.rio.to_raster(base_dir.joinpath("ndvi_clipped.tif"))

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.