<jupyter_start><jupyter_text>Workflow part 1___In the first part of workflow we will *clip rasters, normalize, create visible vegetation indexes*, and *resample images*. Along the process, I decied to save each image before some process, as clipping, for example. Probably wasn't the best decision, but since I worked with drone images, I needed to avoid loading all images on memory, to avoid MemoryErrors. Keep in mind that most of classes and functions cited on workflows are in **module_func_rast**.-- Developed for: Laiza Cavalcante de A. S.-- Period: Msc. Engineering Agriculture<jupyter_code>import os
import re
import rasterio
import numpy as np
from tqdm import tqdm
import seaborn as sns
import geopandas as gp
from matplotlib import pyplot as plt
from module_func_rast import *<jupyter_output><empty_output><jupyter_text>1st Step: Clipping UAV rasters by shapefile<jupyter_code># Raster and shp directories
rast_dir = "/path/file.tif"
shp_dir = "/path/file.shp"

# Creating a list of containing files in a given directory
list_rasters = [ os.path.join(rast_dir, raster) for raster in os.listdir(rast_dir) if raster.endswith('.tif') ]

# Iterating over and clipping
for rast in tqdm(list_rasters):
    rast_class = Raster_operations(rast, shp_dir)

    # Clipping
    new_rast = rast_class.clip_raster_by_shp()<jupyter_output><empty_output><jupyter_text>2nd Step: Normalize UAV rasters<jupyter_code># Ignore division by 0 - avoid errors
np.seterr(divide='ignore', invalid='ignore')

# Directories for rasters, shapes and output
rast_dir = '/path_rasters'
out_directory = '/path_rasters'

# Creating a list of containing files in a given directory
list_rasters = [ os.path.join(rast_dir, raster) for raster in os.listdir(rast_dir) if raster.endswith('.tif') ]

# Iterating over and clipping
for rast in tqdm(list_rasters):
    rast_class = UAV_funcs(rast)

    # Normalizing
    new_rast = rast_class.band_normalized_t1()<jupyter_output><empty_output><jupyter_text>3rd Step: Vegetation Index from UAV normalized rasters<jupyter_code>def vegetation_index(rastB, rastG, rastR, iv_type):
    '''
    ### Visible Vegetation Index creation
        Vegetation Indexes: Excess Green (EXG), Green Leaf Index (GLI),
        Green Red Veg. Index (GRVI), Modified GRVI (MGRVI),
        Red Green Blue Veg. Index (RGBVI), Visible Atmospherically Resistant Index (VARI)
        Input: Blue, Green, Red rasters and each iv must be created
    '''

    def correct_vi(array_):
        array_[array_<-1] = -1
        array_[array_>1] = 1
        return array_

    with rasterio.Env(num_threads='all_cpus'):
        with rasterio.open(rastB, num_threads='all_cpus') as b:
            with rasterio.open(rastG, num_threads='all_cpus') as g:
                with rasterio.open(rastR, num_threads='all_cpus') as r:
                    info = b.meta.copy()
                    info.update({'nodata': np.nan, 'count': 1, 'compress': 'lzw', 'dtype': 'float32'})

                    # Reading raster by tiles (raster windows)
                    tiles = r.block_windows(1)

                    out_exg = img[2][0][:-4] + '_EXG.tif'
                    out_gli = img[2][0][:-4] + '_GLI.tif'
                    out_grvi = img[2][0][:-4] + '_GRVI.tif'
                    out_mgrvi = img[2][0][:-4] + '_MGBVI.tif'
                    out_rgbvi = img[2][0][:-4] + '_RGBVI.tif'
                    out_vari = img[2][0][:-4] + '_VARI.tif'

                    if 'EXG' == iv_type:
                        print('Creating EXG index')
                        with rasterio.open(out_exg, 'w', **info) as dst:
                            for idx, window in tiles:   
                                b_band = b.read(1, window=window, masked=True)
                                g_band = g.read(1, window=window, masked=True)
                                r_band = r.read(1, window=window, masked=True)         
                            
                                iv_calc = ( (2 * g_band) - r_band - b_band )                          
                                dst.write_band(1, iv_calc, window=window)

                    if 'GLI' == iv_type:
                        print('Creating GLI index')
                        with rasterio.open(out_gli, 'w', **info) as dst:
                                for idx, window in tiles:   
                                    b_band = b.read(1, window=window, masked=True)
                                    g_band = g.read(1, window=window, masked=True)
                                    r_band = r.read(1, window=window, masked=True) 

                                    iv_calc = ( ((2 * g_band) - r_band - b_band) /((2 * g_band) + r_band + b_band) )    
                                    iv_calc = correct_vi(iv_calc)                     
                                    dst.write_band(1, iv_calc, window=window)

                    if 'GRVI' == iv_type:
                        print('Creating GRVI index')
                        with rasterio.open(out_grvi, 'w', **info) as dst:
                            for idx, window in tiles:   
                                b_band = b.read(1, window=window, masked=True)
                                g_band = g.read(1, window=window, masked=True)
                                r_band = r.read(1, window=window, masked=True)         
    
                                iv_calc = ((g_band - r_band)/(g_band + r_band))   
                                iv_calc = correct_vi(iv_calc)                     
                                dst.write_band(1, iv_calc, window=window)

                    if 'MGRVI' == iv_type:
                        print('Creating MGRVI index')
                        with rasterio.open(out_mgrvi, 'w', **info) as dst:
                            for idx, window in tiles:   
                                b_band = b.read(1, window=window, masked=True)
                                g_band = g.read(1, window=window, masked=True)
                                r_band = r.read(1, window=window, masked=True)         
                            
                                iv_calc = ( ((g_band**2) - (r_band**2))/((g_band**2) + (r_band**2)) )                           
                                iv_calc = correct_vi(iv_calc)     
                                dst.write_band(1, iv_calc, window=window)

                    if 'RGBVI' == iv_type:
                        print('Creating RGBVI index')
                        with rasterio.open(out_rgbvi, 'w', **info) as dst:
                            for idx, window in tiles:   
                                b_band = b.read(1, window=window, masked=True)
                                g_band = g.read(1, window=window, masked=True)
                                r_band = r.read(1, window=window, masked=True)         
                            
                                iv_calc = ((g_band**2)-(b_band * r_band))/((g_band**2)+(b_band * r_band))                          
                                iv_calc = correct_vi(iv_calc)     
                                dst.write_band(1, iv_calc, window=window)

                    if 'VARI' == iv_type:
                        print('Creating VARI index')
                        with rasterio.open(out_vari, 'w', **info) as dst:
                            for idx, window in tiles:   
                                b_band = b.read(1, window=window, masked=True)
                                g_band = g.read(1, window=window, masked=True)
                                r_band = r.read(1, window=window, masked=True)         
                            
                                iv_calc = ((g_band - r_band)/(g_band + r_band - b_band))                        
                                iv_calc = correct_vi(iv_calc)     
                                dst.write_band(1, iv_calc, window=window)
    
    return iv_calc

rast_dir = '/path_rasters'
np.seterr(divide='ignore', invalid='ignore')

lst_rast = [ os.path.join(rast_dir, raster) for raster in os.listdir(rast_dir) if raster.endswith('.tif') ]
dates = set([re.search(r"(\d{8}|\d{4}-\d{2}-\d{2}|\d{2}-\d{2}-\d{4}|\d{4}_\d{2}_\d{2})", elm).group() for elm in lst_rast])

rast = []
for date in dates:
    bands, elms = [], []
    for elm in lst_rast:
        if date in elm:
            bands.append(re.search(r"([R|G|B]+)", elm).group())
            elms.append(elm)
    rast.append( (date, bands, elms) )

ivs = ['GRVI', 'GLI', 'VARI', 'RGBVI', 'EXG']

for img in tqdm(rast):
    for iv in tqdm(ivs):
        result = vegetation_index(rastB = img[2][0], 
                        rastG = img[2][1], 
                        rastR = img[2][2],
                        iv_type = iv)<jupyter_output><empty_output><jupyter_text>4th Step: Resample images<jupyter_code>rast_dir = '/path_rasters'
lst_rast = [ os.path.join(rast_dir, raster) for raster in os.listdir(rast_dir) if raster.endswith('.tif') ]

Sentinel = r'/path_raster_base.tif'

# Iterating over and clipping
for rast in tqdm(lst_rast):
    rst_class = Raster_resample(raster_base=Sentinel, 
                                raster_dir=rast, 
                                scale = 132)
    rst_resample = rst_class.resample_by_raster()<jupyter_output><empty_output><jupyter_text>5th Step:Looking for Outliers and Cheking normality<jupyter_code># Testing outlier removal
path = r'/path_rasters/'

# Looking for outliers in array
outlier = Outliers_check(path_img = path)
outliers_list = outlier.find_outlier()
print(outliers_list)

# Removing outliers
new_raster = Outliers_check(path_img = path, outliers=outliers_list).remove_outlier()

# Checking normality
normality = Outliers_check(path_img = path)
rst = normality.normality_check()<jupyter_output><empty_output><jupyter_text>6th Step: Plotting Raster<jupyter_code># Plotting raster with outlier and without
with rasterio.open(path) as tif:
    org_raster = tif.read(1)
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8,4))
fig = ax1.imshow(org_raster, cmap='hot', vmin=-0.2, vmax=0.2)
fig = ax2.imshow(new_raster, cmap='hot', vmin=-0.2, vmax=0.2)

ax1.set_title('Original')
ax2.set_title('Without outlier')
plt.colorbar(fig)

# Plotting raster histogram with outlier and without
sns.set_style("dark")
sns.despine()

array_org = org_raster[np.logical_not(np.isnan(org_raster))].flatten()
array_mod = new_raster[np.logical_not(np.isnan(new_raster))].flatten()

fig2, axes = plt.subplots(nrows=1, ncols=2, figsize=(8,4), sharey=True)
for ax, array in zip(axes, [array_org, array_mod]):
    sns.distplot(array_org, ax=ax, 
            hist_kws={'alpha': 0.4, 'color':'purple'},
            kde_kws={'color': 'purple'})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11254199328107503, 'fraction_numerical': 0.010544979469951475, 'mean_word_length': 1.8411983032873807, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '310969', 'license': 'MIT', 'n_tokens_mistral': 3573, 'n_tokens_neox': 3329, 'n_words': 766, 'path': 'notebooks/01_Workflow.ipynb', 'repo_name': 'cavalcante-l/UAV', 'size': '10716'}