Composites

import glob, os
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
import rasterio.plot
import datetime
from pathlib import Path
from IPython.display import display

print('All libraries successfully imported!')

All libraries successfully imported!

Set directory

computer_path = '/export/miro/ndeffense/LBRAT2104/'
grp_letter    = 'X'

data_path = f'{computer_path}data/'                     # Directory with data shared by the assistant
work_path = f'{computer_path}GROUP_{grp_letter}/WORK/'  # Directory for all work files

# Image timeserie directory
ts_path = f'{work_path}NDVI/'
raster_file_list = sorted(glob.glob(f'{ts_path}*.tif'))

# Output directory
composite_path = f'{work_path}COMPOSITE/'

Path(composite_path).mkdir(parents=True, exist_ok=True)

Set parameters

Compositing method :

• mean
• median
• maximum
• minimum

method = 'minimum'

Date

The date must be in each image filename. begin_date is the first position of the date

E.g. : "SITE320181108_NDVI.tif" --> begin_date = 8

format     = '%Y%m%d' # format : YYYYMMDD
begin_date = 8

Set filenames

composite_tif = f'{composite_path}composite_{method}.tif'

Get dataframe with all images

• date of acquisition
• file path

dict_list = []

for im in raster_file_list:
            
        date_str = os.path.basename(im)[begin_date-1:begin_date-1+8]

        date_im = datetime.datetime.strptime(date_str, format).date()

        dict_list.append({'date': date_im,
                          'im_path': im})

im_date_path_df = pd.DataFrame.from_dict(dict_list).sort_values('date')

display(im_date_path_df)

	date	im_path
0	2020-01-16	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
1	2020-02-12	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
2	2020-03-16	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
3	2020-04-17	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
4	2020-05-20	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
5	2020-06-21	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
6	2020-07-19	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
7	2020-08-13	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
8	2020-09-14	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
9	2020-10-19	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
10	2020-11-18	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
11	2020-12-18	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...

Select images inside a period

start_date = datetime.datetime.strptime('20200101', format).date()
end_date   = datetime.datetime.strptime('20200601', format).date()

mask = (im_date_path_df['date'] > start_date) & (im_date_path_df['date'] <= end_date)

im_date_path_df = im_date_path_df.loc[mask]

display(im_date_path_df)

	date	im_path
0	2020-01-16	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
1	2020-02-12	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
2	2020-03-16	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
3	2020-04-17	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...
4	2020-05-20	/export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/...

list_src_arr = []

for i, row in im_date_path_df.iterrows():
    im_file    = row['im_path']
    src = rasterio.open(im_file, "r")
    im_arr = src.read(1)
    list_src_arr.append(im_arr)
    src.close()
    
print(f'Shape of features : {im_arr.shape}')
print(f'Number of features : {len(list_src_arr)}')

Shape of features : (570, 986)
Number of features : 5

im_arr_stack = np.dstack(list_src_arr)

print(im_arr_stack.shape)
print(f'There are {im_arr_stack.shape[2]} features')
print(f'The features type is : {im_arr_stack.dtype}')

(570, 986, 5)
There are 5 features
The features type is : float64

if method == 'median':
    compo_arr = np.nanmedian(im_arr_stack, axis=2)

elif method == 'mean':
    compo_arr = np.nanmean(im_arr_stack, axis=2)

elif method == 'std':
    compo_arr = np.nanstd(im_arr_stack, axis=2)

elif method == 'minimum':
    compo_arr = np.nanmin(im_arr_stack, axis=2)

elif method == 'maximum':
    compo_arr = np.nanmax(im_arr_stack, axis=2)

else:
    print(f'Method not available !')

print(compo_arr)

[[0.42912925 0.43237788 0.37547408 ... 0.5170441  0.49381541 0.34527089]
 [0.43097239 0.43376623 0.42655262 ... 0.49025202 0.32620818 0.20146688]
 [0.43283582 0.44285215 0.44069015 ... 0.33983287 0.18755803 0.14951095]
 ...
 [0.59484778 0.55236576 0.55077539 ... 0.32256267 0.24806746 0.31261561]
 [0.60609911 0.56308411 0.55648536 ... 0.32063146 0.24914676 0.25447244]
 [0.66213314 0.596366   0.59649123 ... 0.3208756  0.23613963 0.26780627]]

/export/miro/ndeffense/.local/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: All-NaN slice encountered
  # This is added back by InteractiveShellApp.init_path()

img_temp_tif = raster_file_list[0]

with rasterio.open(img_temp_tif) as src:
    profile = src.profile

print(profile)

with rasterio.open(composite_tif, "w", **profile) as dst:
    dst.write(compo_arr, 1)

{'driver': 'GTiff', 'dtype': 'float64', 'nodata': -10000.0, 'width': 986, 'height': 570, 'count': 1, 'crs': CRS.from_epsg(32631), 'transform': Affine(10.0, 0.0, 627260.0,
       0.0, -10.0, 5596180.0), 'tiled': False, 'compress': 'lzw', 'interleave': 'band'}