Skip to content
Extraits de code Groupes Projets
Valider a1e6579e rédigé par Dries De Bièvre's avatar Dries De Bièvre
Parcourir les fichiers

preprocessing steps in seperate functions

parent 9ef200c5
Aucune branche associée trouvée
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
import glob, os
import rasterio
import numpy as np
from rasterio.enums import Resampling
def resample_SCL(list_L2A, resampled_path, upscale_factor):
for L2A_safe in list_L2A:
print(f'Start resampling SCL of L2A image: {L2A_safe}')
im_file_20m = glob.glob(f'{L2A_safe}/GRANULE/*/IMG_DATA/R20m/*_SCL_20m.jp2')[0]
im_file_20m = im_file_20m.replace('\\','/')
im_file_10m = f'{resampled_path}{os.path.basename(im_file_20m).replace("20m","10m").replace("jp2","tif")}'
if not os.path.isfile(im_file_10m):
# Open file
src = rasterio.open(im_file_20m, "r", driver='JP2OpenJPEG')
# Resample data to target shape
resampled_data = src.read(out_shape = (src.count,
int(src.height * upscale_factor),
int(src.width * upscale_factor)
),
resampling = Resampling.nearest)
# Scale image transform
new_transform = src.transform * src.transform.scale(
(src.width / resampled_data.shape[-1]),
(src.height / resampled_data.shape[-2])
)
# Update metadata
profile = src.profile
profile.update(driver = 'GTiff',
width = src.width*upscale_factor,
height = src.height*upscale_factor,
transform = new_transform)
# Write resampled image
dst = rasterio.open(im_file_10m, "w", **profile)
dst.write(resampled_data)
# Close rasterio objects
src.close()
dst.close()
print(f'--> A new resampled raster file is created : {im_file_10m}')
print('--> All SCL are resampled !')
def resample_bands(list_L2A, resampled_path, upscale_factor, bands_20m_list, resampling_method=Resampling.bilinear):
for L2A_safe in list_L2A:
print(f'Start resampling L2A image: {L2A_safe}')
for band in bands_20m_list:
print("Start resampling band : ", band)
im_file_20m = glob.glob(f'{L2A_safe}/GRANULE/*/IMG_DATA/R20m/*{band}_20m.jp2')[0]
im_file_20m = im_file_20m.replace('\\','/')
im_file_10m = f'{resampled_path}{os.path.basename(im_file_20m).replace("20m","10m").replace("jp2","tif")}'
if not os.path.isfile(im_file_10m):
# Open file
src = rasterio.open(im_file_20m, "r", driver='JP2OpenJPEG')
# Resample data to target shape
resampled_data = src.read(out_shape = (src.count,
int(src.height * upscale_factor),
int(src.width * upscale_factor)
),
resampling = resampling_method)
# Scale image transform
new_transform = src.transform * src.transform.scale(
(src.width / resampled_data.shape[-1]),
(src.height / resampled_data.shape[-2])
)
# Update metadata
profile = src.profile
profile.update(driver = 'GTiff',
width = src.width*upscale_factor,
height = src.height*upscale_factor,
transform = new_transform)
# Write resampled image
dst = rasterio.open(im_file_10m, "w", **profile)
dst.write(resampled_data)
# Close rasterio objects
src.close()
dst.close()
print(f'--> A new resampled raster file is created : {im_file_10m}')
print('--> All reflectances are resampled !')
def clip_imgs(list_im_to_clip, clipped_path, roi_gdf):
for im_file in list_im_to_clip:
im_file_roi = f'{clipped_path}{os.path.basename(im_file)[:-4]}_ROI.tif'
if not os.path.isfile(im_file_roi):
# Open file
src = rasterio.open(im_file, "r")
# Clip the raster to the extent of the shape
out_image, out_transform = rasterio.mask.mask(src,
roi_gdf.geometry,
all_touched=True,
crop=True)
# Update metadata
out_meta = src.meta
out_meta.update(driver='GTiff',
width = out_image.shape[2],
height = out_image.shape[1],
transform = out_transform)
# Write clipped image
dst = rasterio.open(im_file_roi, "w", **out_meta)
dst.write(out_image)
# Close rasterio objects
src.close()
dst.close()
print(f'A new raster file is created : {im_file_roi}')
print('--> All images are clipped !')
def apply_SCL(list_im_ROI, clipped_path, masked_path, nodata_val):
for im_file in list_im_ROI:
# Get date of image
date = os.path.basename(im_file)[7:7+15]
# Get tile of image
tile = os.path.basename(im_file)[0:6]
# Find SCL corresponding to the given reflectances image (same date and same tile)
scl_file = glob.glob(f'{clipped_path}*{tile}*{date}*SCL_10m_ROI.tif')[0]
scl_file = scl_file.replace('\\','/')
im_file_scl = f'{masked_path}{os.path.basename(im_file)[:-4]}_SCL.tif'
if not os.path.isfile(im_file_scl):
# Open SCL and change invalid pixels categories by NaN
src = rasterio.open(scl_file, "r")
# Read file as numpy array
SCL = src.read(1)
src.close()
#print('Scene Classification map')
#show(SCL, cmap='Set3')
SCL = SCL.astype(float)
SCL[SCL == 0] = np.nan # No data
SCL[SCL == 1] = np.nan # Saturated or defective
SCL[SCL == 2] = np.nan # Dark area pixels
SCL[SCL == 3] = np.nan # Cloud shadows
SCL[SCL == 4] = 1 # Vegetation
SCL[SCL == 5] = 1 # Not vegetated
SCL[SCL == 6] = 1 # Water
SCL[SCL == 7] = 1 # Unclassified
SCL[SCL == 8] = np.nan # Cloud medium probability
SCL[SCL == 9] = np.nan # Cloud high probability
SCL[SCL == 10] = np.nan # Thin cirrus
SCL[SCL == 11] = np.nan # Snow
# Open file
src = rasterio.open(im_file, "r")
# Read file as numpy array
im = src.read(1)
# Update metadata
profile = src.profile
profile.update(dtype=rasterio.int16, # Set to int16 it is lighter than float
nodata=nodata_val, # Set nodata value in metadata
compress='lzw') # Compression option
# Mask image reflectance with SCL
im_SLC = im * SCL
# Change numpy NaN by nodata_val (e.g. -10000)
im_SLC[np.isnan(im_SLC)] = nodata_val
# Change the array's type : from float to integer 16
im_SLC = im_SLC.astype(np.int16)
# Write image
dst = rasterio.open(im_file_scl, 'w', **profile)
dst.write(im_SLC, 1)
# Close rasterio objects
src.close()
dst.close()
print(f'A new raster file is created : {im_file_scl}')
print('--> SCL is applied on all images !')
\ No newline at end of file
Impossible d'afficher diff de source : il est trop volumineux. Options pour résoudre ce problème : voir le blob.
0% Chargement en cours ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter