#python code to calculate attributable flood depth change during Hurricane Harvey
# by Michael F Wehner 04.26.21
# Lawrence Berkeley National Laboratory
#
# Usage: python ./make_Harvey_flood_difference.py reference_file.tif climate_change_file.tif difference_file.tif
# Where:
# reference_file.tif is usually the actual flood geotiff file (i.e. 0percent_human_effect.tif)
# climate_change_file.tif is a selected amount of human influence geotiff file (i.e. 38percent_human_effect.tif)
# difference_file.tif is the output geotiff file containing the attributable change in flood depth (meters)


import gdal,sys
def getRasterBand(fn,band=1,access=0):
 ds=openRaster(fn,access)
 band=ds.GetRasterBand(1).ReadAsArray()
 return band

def openRaster(fn,access=0):
 ds=gdal.Open(fn,access)
 return ds

def createRasterFromCopy(fn,ds,data,driverFmt="GTiff"):
 driver=gdal.GetDriverByName(driverFmt)
 outds=driver.CreateCopy(fn,ds,strict=0)
 outds.GetRasterBand(1).WriteArray(data)
 ds=None
 outds=None

f1=sys.argv[1]
f8=sys.argv[2]
fout=sys.argv[3]

var1=getRasterBand(f1)
var8=getRasterBand(f8)
diff18=var1-var8

driver=gdal.GetDriverByName("GTiff")
outds=driver.CreateCopy(fout,gdal.Open(f1),strict=0)
outds.GetRasterBand(1).WriteArray(diff18)
outds.None