ROL

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

#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]

#get the mean values for W-E.
WE_mean=[]
for j in range (len(data)):
    temp_mean=np.mean(data[j])
    WE_mean.append(temp_mean)
    
thsd_mean=np.mean(WE_mean)
thsd_std=np.mean(WE_mean)
thsd=thsd_mean+thsd_std

for j in range (len(WE_mean)):
    if WE_mean[j]>thsd:
        print(j)