#!/bin/bash -x

#input is yearly GWP file
input=$1
indir=$2
outdir=$3

export GDAL_CACHEMAX=6000

date=$(basename $input .tif | cut -d'_' -f2)

yr=${date:0:4}

dir_to_search=${indir}/ocean_mask/
cp ${dir_to_search}/ocean_mask_global_GWP_res.tif .
ls .

if [ ! -f ${outdir}/GWP.OSWF.ANNUAL.${yr}.v1.tif ];
then
	mkdir -p output/
	
	#reproject to WGS84
	gdalwarp -t_srs EPSG:4326 $input swat_${yr}_global.tif -te -180 -90 180 90 -ts 172800 86400 -co compress=deflate -r near -dstnodata 65535 -ot UInt16
	#change nodata value
	gdal_translate swat_${yr}_global.tif swat_${yr}_global_1.tif -a_nodata 65534 -ot UInt16 -co compress=deflate	
	
	#account for leap years
	if [[ $yr == "2000" || $yr == "2004" || $yr == "2008" || $yr == "2012" || $yr == "2016" || $yr == "2020" || $yr == "2024" ]];
	then
		#change pixel value 65535 to 0
		gdal_calc.py --outfile=swat_${yr}_global_2.tif -A swat_${yr}_global_1.tif --calc="0*(A==65535)+A*logical_and(A>0, A<65535)" --co=COMPRESS=LZW --type=UInt16
		#add ocean mask to GWP file
		gdal_calc.py --outfile=swat_masked_${yr}.tif -A swat_${yr}_global_2.tif -B ocean_mask_global_GWP_res.tif --calc="A*B+366*(B==0)" --co=COMPRESS=LZW --type=UInt16
	else
		gdal_calc.py --outfile=swat_${yr}_global_2.tif -A swat_${yr}_global_1.tif --calc="0*(A==65535)+A*logical_and(A>0, A<65535)" --co=COMPRESS=LZW --type=UInt16
		gdal_calc.py --outfile=swat_masked_${yr}.tif -A swat_${yr}_global_2.tif -B ocean_mask_global_GWP_res.tif --calc="A*B+365*(B==0)" --co=COMPRESS=LZW --type=UInt16
	fi
	
	#change tif to cog file
	gdal_translate -of COG -co COMPRESS=ZSTD -co PREDICTOR=2 -co RESAMPLING=NEAREST swat_masked_${yr}.tif GWP.OSWF.ANNUAL.${yr}.v1.tif
	
	mv GWP.OSWF.ANNUAL.${yr}.v1.tif output/
else
	echo "${outdir}/GWP.OSWF.ANNUAL.${yr}.v1.tif already exists"
fi

exit 0