10.fragmentation_indicatior_cal.py

# -*- coding: utf-8 -*-
"""
Created on Sat Nov 11 12:39:02 2023

@author: no67wuwu
"""

import os
import numpy as np
import pandas as pd
from glob import glob

import MacPyver as mp

from multiprocessing import Pool

# fragmentation function from Aurora
def frag_ind(s):
    area = ((sur_radius*2)+1)**2
    a = map(lambda x: x**2, s)
    a = sum(a)*1.0/area
    return a


# path to files:
tile_path = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\result_grass\frag_unique_masled_tr1000"


# path to textfile / split the file / result from the grass report function
report_path = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\result_grass\tables\report_unique_areas.file"
df_report = pd.read_csv(report_path, sep='|', skiprows=15 ,skipfooter=3, names=['a', 'frag_id', 'b','area','c'], thousands=',', engine='python')
df_report = df_report.loc[:, ['frag_id', 'area']]

outpath = r'S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\result\frag_ind'
log = open(outpath+'\\frag_ind.log', 'w')



def get_tiles(tile):
    ''' get surrounding tiles '''
    name = tile.split(os.sep)[-1]
    y_top = int(name.split('-')[1])-1
    y_bot = y_top + 2
    x_left = int(name.split('-')[2].split('.')[0])-1
    x_right = x_left + 2
    
    #all possible surrounding tiles
    sur_tiles = ['fr_cl_tile-{}-{}.tif'.format(('000'+str(y))[-3:], ('000'+str(x))[-3:]) 
            for y in range(y_top, y_bot+1)  
                for x in range(x_left, x_right+1) ]
     
    #position of the tiles               
    pos_dic_tiles = {}
    pos = ['0,0','0,1','0,2',
           '1,0','1,1','1,2',
           '2,0','2,1','2,2']
           
    pos_count = 0
    for st in sur_tiles:
        pos_dic_tiles[st] = pos[pos_count]
        pos_count +=1
    
    
    #get all existing tiles
    sur_tile_paths = [x for x in tiles for y in sur_tiles if y in x]
    
    # create empty array
    array = np.zeros([300,300])
    array.fill(-1)
    
    # fill array with data
    for stp in sur_tile_paths:
        # get position in merge array
        stp_name = stp.split(os.sep)[-1]
        sub_rast = mp.raster.tiff.read_tif(stp)
        
        # Check the shape of the sub_rast
        if sub_rast.shape != (100, 100):
            print(f"Unexpected shape for {stp_name}: {sub_rast.shape}")
            continue  # Skip this tile or handle the error as needed

        # fill position based on the position dic
        ymin = int(pos_dic_tiles[stp_name].split(',')[0]) * 100
        ymax = ymin +100
        xmin = int(pos_dic_tiles[stp_name].split(',')[1]) * 100
        xmax = xmin + 100
        array[ymin:ymax, xmin:xmax] = sub_rast
    
    return array
    
    # # create empty array
    # array = np.zeros([300,300])
    # array.fill(-1)
    # #fill array with data
    
    # for stp in sur_tile_paths:
    #     #get position in merge array
    #     stp_name = stp.split(os.sep)[-1]
    #     sub_rast = mp.raster.tiff.read_tif(stp)
    #     # fill position based on the position dic
    #     ymin = int(pos_dic_tiles[stp_name].split(',')[0]) * 100
    #     ymax = ymin +100
    #     xmin = int(pos_dic_tiles[stp_name].split(',')[1]) * 100
    #     xmax = xmin + 100
    #     array[ymin:ymax, xmin:xmax] = sub_rast
    
    # return array
    
    
    
def bbox(x,y, sur_radius):
    '''get the bounding box'''
    xmin = x - sur_radius
    xmax = x + sur_radius + 1
    ymin = y - sur_radius
    ymax = y + sur_radius + 1
    return xmin, xmax, ymin, ymax

def main(tile):
    
    # get surrownding tiles
    raster = get_tiles(tile)
    orig_raster = mp.raster.tiff.read_tif(tile)
    out_count_area = np.zeros(orig_raster.shape)
    out_orig_area = np.zeros(orig_raster.shape)
    
    frag_ids = np.unique(raster[:,:]) 
    #frag_ids = np.unique(raster[95:205, 95:205]) 
    frag_ids = np.delete(frag_ids, np.where(frag_ids==0))
    df_rep_subset = df_report.loc[df_report['frag_id'].isin(frag_ids)]
    
    
    for x in range(100,200):
        for y in range(100, 200):
            if raster[x,y]!=0:
                # get surrounding cells
                xmin, xmax, ymin, ymax = bbox(x, y, sur_radius)
                s_cells = raster[ xmin : xmax, ymin : ymax]
                window_ids, window_counts = np.unique(s_cells, return_counts=True)
                if 0 in window_ids:
                    window_ids = window_ids[1:]
                    window_counts = window_counts[1:]
                #calc fragmentation based on pixel count in moving window
                frag_indi_counts = frag_ind(window_counts)
                
                #calc fragmentation based on orig areas
                df_window = df_rep_subset.loc[df_rep_subset['frag_id'].isin(window_ids)]
                df_window_array = np.array(df_window['area']/(1000*1000)).astype(int)
                frag_indi_orig = frag_ind(df_window_array)
                
                
                out_count_area[x-100,y-100] = int(frag_indi_counts + .5)
                out_orig_area[x-100,y-100] = int(frag_indi_orig + .5)

    
    # considering just the pixel in the window    
    outraster_counts = r'S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\result\window_count3' + os.sep + tile.split(os.sep)[-1]
    mp.raster.tiff.write_tif(tile,outraster_counts, out_count_area, 1)
    # considering the original area from the id    
    outraster_orig_area = r'S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\result\orig_area3' + os.sep + tile.split(os.sep)[-1]
    mp.raster.tiff.write_tif(tile,outraster_orig_area, out_orig_area, 1)
    return 'tile: {} written'.format(tile.split(os.sep))


tiles = glob(tile_path + '\*.tif')
sur_radius = 5 #radius of cells for surrounding the center_pixel

if __name__=='__main__':
    pool = Pool(20)
    for bar in pool.imap_unordered(main, tiles):
        log.write(bar)
        log.flush()

    log.close()