AIL

AIL

#AIL
#Version 1.0
#Xinyuan Wei 03/19/2019
#Edge Identification from the Simulated TSLF Map
'''
This program can identify the edge from a simulated TSLF map.
The identification algorithm is Agglomerating Partial Inner Lines (APIL).
'''
import gdal
import osr
import os
import ogr
import math
import numpy as np
import scipy.stats
from gdalconst import *
gdal.AllRegister()

#size of the center (cells)
center_width=2100 

#get the work directory
print('[The work directory:]',os.getcwd())

#read the raster data
dataname = 'Grid_Burned_Times_1200.img'
rasterdata = gdal.Open(dataname,GA_ReadOnly)
if rasterdata is None:
    print ('Could not open ' + dataname)
    sys.exit(1)

#get the raster data information
cols = rasterdata.RasterXSize
rows = rasterdata.RasterYSize
bands = rasterdata.RasterCount
driver = rasterdata.GetDriver().LongName
print('')
print('[Description of the study area map:]')
print('Size of the map: '+str(cols)+'*'+str(rows)+' ha'+'.')
print('Number of bands: '+str(bands)+'.')
print('Format of this image: '+str(driver)+'.')

geotransform = rasterdata.GetGeoTransform()
print('')
print('xOrigin,pixelWidth,xrotation,yOrigin,xrotation,pixelHeight')
print(geotransform)

#get the band information
band=rasterdata.GetRasterBand(1)
bandtype=gdal.GetDataTypeName(band.DataType)
print('')
print('[Band type is]'+bandtype+'.')

#convert the image data to an array
data=band.ReadAsArray(0,0,cols,rows)
print('[Burned times of each pixel (year)]:')
print(data)
print('')
print('[Shape of the data]:',np.shape(data))

for i in range (len(data)):
    for j in range (len(data[0])):
        data[i][j]=1/data[i][j]

#calculated the mean of the center 2100*2100
mean=np.mean(data[400:2499,400:2499])
p_std=((np.std(data[400:2499,400:2499]))*
       (math.sqrt(len(data[400:2499,400:2499])*len(data[400:2499,400:2499][0]))))
print('')
print('[The mean of the center size as 2100*2100]:',mean)
print('[The population standard deviation of the center size as 2100*2100]:',p_std)
print((np.std(data))*(math.sqrt(len(data)*len(data[0]))))
#get the center area
center=data[400:2500,400:2500]
print('')
#print("The center is")
#print(center)
print('[Shape of the data]:',np.shape(center))

#agglomerative calculation
a=round((2900-center_width)/2)
b=round((2900-center_width)/2+center_width)
c=round((2900-center_width)/2)
d=round((2900-center_width)/2+center_width)
#print(a,b,c,d)

std_array=[]

for i in range(10000):
    std1=np.std(data[a-1:b,c:d])*(math.sqrt(len(data[a-1:b,c:d])*len(data[a-1:b,c:d][0]))) #top
    std2=np.std(data[a:b+1,c:d])*(math.sqrt(len(data[a:b+1,c:d])*len(data[a:b+1,c:d][0]))) #botom
    std3=np.std(data[a:b,c-1:d])*(math.sqrt(len(data[a:b,c-1:d])*len(data[a:b,c-1:d][0]))) #left
    std4=np.std(data[a:b,c:d+1])*(math.sqrt(len(data[a:b,c:d+1])*len(data[a:b,c:d+1][0]))) #right
    std=(std1,std2,std3,std4)
    
    order_n=std.index(min(std))
    
    if order_n==0:
        a=a-1
    if order_n==1:
        b=b+1
    if order_n==2:
        c=c-1
    if order_n==3:
        d=d+1
        
    if a==0:
        break
    if b==2899:
        break
    if c==0:
        break
    if d==2899:
        break
    std_array.append(min(std))
    #print(a,b,c,d)
    print(min(std))

#the array without edge area   
array_new=(data[a:b,c:d])
    
print('')
print(array_new)
print('[Shape of the array:]')
np.savetxt('AIL_Map.csv',array_new,delimiter=',')
print(np.shape(array_new))
np.savetxt('AIL_SD.csv',std_array,delimiter=',')