#Extract Raster Change Areas
#Chance Morrison, Jun 2026

import arcpy, os

#Set workspace and get parameters

#ws = os.getcwd()
#arcpy.env.scratchWorkspace = 'memory'
#mem_pth = 'memory/ext_ras_chng_area'

gdb = arcpy.GetParameterAsText(0)
ras1 = arcpy.GetParameterAsText(1)
ras0 = arcpy.GetParameterAsText(2)
ply_out = arcpy.GetParameterAsText(3)
con_where = arcpy.GetParameter(4)
fs_size = arcpy.GetParameter(5)
chng_small = arcpy.GetParameter(6)
chng_large = arcpy.GetParameter(7)
ras_save = arcpy.GetParameter(8)
ras_out = arcpy.GetParameterAsText(9)

arcpy.env.workspace = gdb

#Compare raster details before processing

ras0_sr = arcpy.Describe(ras0).spatialReference
ras1_sr = arcpy.Describe(ras1).spatialReference

arcpy.AddMessage(f'{ras0_sr.name}')
arcpy.AddMessage(f'{ras1_sr.name}')

if ras0_sr.name == ras1_sr.name:
    arcpy.AddMessage('The spatial references match')
else:
    arcpy.AddMessage('The spatial references don\'t match')

ras0_desc = arcpy.Describe(ras0)
ras1_desc = arcpy.Describe(ras1)

arcpy.AddMessage(f'{ras0_desc.bandCount}')
arcpy.AddMessage(f'{ras1_desc.bandCount}')

if ras0_desc.bandCount == ras1_desc.bandCount:
    arcpy.AddMessage('The band counts match')
else:
    arcpy.AddMessage('The band counts don\'t match')

arcpy.AddMessage(f'{ras0_desc.meanCellHeight}')
arcpy.AddMessage(f'{ras1_desc.meanCellHeight}')

if ras0_desc.meanCellHeight == ras1_desc.meanCellHeight:
    arcpy.AddMessage('The cell heights match')
else:
    arcpy.AddMessage('The cell heights don\'t match')

arcpy.AddMessage(f'{ras0_desc.meanCellWidth}')
arcpy.AddMessage(f'{ras1_desc.meanCellWidth}')

if ras0_desc.meanCellWidth == ras1_desc.meanCellWidth:
    arcpy.AddMessage('The cell widths match')
else:
    arcpy.AddMessage('The cell widths don\'t match')

arcpy.AddMessage(f'{ras0_desc.extent}')
arcpy.AddMessage(f'{ras1_desc.extent}')

if ras0_desc.meanCellWidth == ras1_desc.meanCellWidth:
    arcpy.AddMessage('The extents match')
else:
    arcpy.AddMessage('The extents don\'t match')

x_diff = abs(ras0_desc.extent.XMin - ras1_desc.extent.XMin)
y_diff = abs(ras0_desc.extent.YMin - ras1_desc.extent.YMin)

arcpy.AddMessage(f'{x_diff}')
arcpy.AddMessage(f'{y_diff}')

arcpy.AddMessage(f'{x_diff % int(ras0_desc.meanCellWidth)}')
arcpy.AddMessage(f'{y_diff % int(ras0_desc.meanCellHeight)}')

if ((x_diff % int(ras0_desc.meanCellWidth)) == 0) and (y_diff % int(ras0_desc.meanCellHeight) == 0):
    arcpy.AddMessage('The grids match')
else:
    arcpy.AddMessage('The grids don\'t match')

if ras_save == True:
    arcpy.AddMessage('Saving each processed raster layer for review...')
else:
    arcpy.AddMessage('Not saving each processed raster layer...')

#Compute relative change raster
arcpy.AddMessage('Calculating relative difference...')
ras_reldif = arcpy.ia.ComputeChangeRaster(ras1, ras0, 'DIFFERENCE')

if ras_save == True:
    ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_1ras_reldif.tif'
    ras_reldif.save(ras_pth)
    arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

ras_in = ras_reldif

#Absolute raster
arcpy.AddMessage('Calculating absolute change...')
ras_abs = arcpy.sa.Abs(ras_in)

if ras_save == True:
    ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_2ras_abs.tif'
    ras_abs.save(ras_pth)
    arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

ras_in = ras_abs

#Con raster for threshold
if con_where != '':
    arcpy.AddMessage(f'Calculating value threshold where "{con_where}"...')
    ras_con = arcpy.sa.Con(ras_in, ras_in, 0, con_where)

    if ras_save == True:
        ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_3ras_con.tif'
        ras_con.save(ras_pth)
        arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

    ras_in = ras_con

#Focal statistics on raster
if fs_size > 0:
    arcpy.AddMessage(f'Reducing noise with a square the size of "{fs_size}"...')
    fs_nbhd = arcpy.sa.NbrRectangle(fs_size, fs_size, 'CELL')
    ras_fs = arcpy.sa.FocalStatistics(ras_in, fs_nbhd, 'MINIMUM')

    if ras_save == True:
        ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_4ras_fs.tif'
        ras_fs.save(ras_pth)
        arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

    ras_in = ras_fs

#Con raster for boolean
arcpy.AddMessage('Calculating boolean change...')
ras_bool = arcpy.sa.Con(ras_in, 1, 0, 'Value > 0')

if ras_save == True:
    ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_5ras_bool.tif'
    ras_bool.save(ras_pth)
    arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

ras_in = ras_bool

#Group regions
arcpy.AddMessage('Grouping change regions...')
ras_reg = arcpy.sa.RegionGroup(ras_in)

if ras_save == True:
    ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_6ras_reg.tif'
    ras_reg.save(ras_pth)
    arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

ras_in = ras_reg

#Con raster for size
#Cannot use SQL query field again since it requires a depenency on a layer,
#and this layer is not created until after the tool is ran...

chng_large_where = ''
chng_small_where = ''
if chng_large > 0:
    chng_large_where = f' And Count <= {chng_large}'
if chng_small > 0:
    chng_small_where = f' And Count >= {chng_small}'
chng_where = f'{chng_large_where}{chng_small_where}'

if chng_where != '':
    ras_chng = arcpy.sa.Con(ras_in, ras_bool, 0, chng_where)

    if ras_save == True:
        ras_pth = ras_out + '\\' + os.path.basename(ply_out) + '_7ras_chng.tif'
        ras_chng.save(ras_pth)
        arcpy.AddMessage(f'"{os.path.basename(ras_pth)}" saved.')

    ras_in = ras_chng

#Raster to polygon
arcpy.AddMessage('Converting raster change groups to polygons...')
arcpy.conversion.RasterToPolygon(ras_in, ply_out, 'NO_SIMPLIFY')

#Remove non-change areas
arcpy.AddMessage('Isolating change polygons...')
with arcpy.da.UpdateCursor(ply_out, ['gridcode'], 'gridcode = 0') as cursor:
    for row in cursor:
        cursor.deleteRow()

#Zonal statistics of ras_reldif within ply_ras
arcpy.AddMessage('Calculating statistics within change polygons...')
arcpy.sa.ZonalStatisticsAsTable(ply_out, 'Id', ras_reldif, 'tbl_sum', 'DATA', 'ALL')
arcpy.management.JoinField(ply_out, 'Id', 'tbl_sum', 'Id')

#Updating polygon layer symbology
#.lyrx file must be in same directory as script file
#arcpy.AddMessage(f'Updating symbology of "{os.path.basename(ply_out)}"')

#lyrx_pth = os.path.join(os.getcwd(), 'raster_change.lyrx')
#aprx = arcpy.mp.ArcGISProject("CURRENT")
#mp = aprx.listMaps()[0]

#lyr = mp.addDataFromPath(ply_out)
#flds = "VALUE_FIELD SUM SUM"
#arcpy.management.ApplySymbologyFromLayer(lyr, lyrx_pth, flds, 'MAINTAIN')
#arcpy.SetParameter(1, lyr)

#Script end
arcpy.AddMessage('Completed script')