Skip to content
Snippets Groups Projects
Select Git revision
  • 0fd9a78835f6807debbf2bee45dd8b7585ee9f2e
  • main default protected
  • develop
  • v0.1
4 results

data_raw.sav

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    10.fragmentation_indicatior_cal.py 5.72 KiB
    # -*- 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()