#!/usr/bin/env perl

########## Check or modify the path to OSGeo4w.bat a
$ogr="C:/Program Files/QGIS 3.34.2/OSGeo4w.bat";
###############################################################

### USER INPUT ################################################
$startyear=2020;
$endyear=2023;
$source="processed_tiles";
###############################################################

########## Check name of the tile list
open (DAT, "tiles.txt") or die "ERROR: no tiles.txt"; @tilelist=<DAT>; close (DAT); foreach (@tilelist) {chomp;}
######################################

$ogr=~s/\\/\//g;

# Setup #
$|=1;
$tileside=4000;
$tilebuffer=2;
$pixelsize=0.00025;
@names=("Loss","TCC","TCH");

for $name (@names){
	for ($year=$startyear; $year<=$endyear; $year++){
		print "aggregating $year $name ...\n";
		$outname="$name\_$year";
		if (!-e "$outname.tif"){
			#single tile
			$tilenum=@tilelist;
			if ($tilenum==1){system "copy $source\\$name\_$year\_$tilelist[0].tif $outname.tif";}
			else {
			#AOI area
			my @harray=(); my @varray=(); 
			for $tile (@tilelist){
			$h=substr $tile,0,3; $h=$h*1; $hsign=substr $tile,3,1; if ($hsign eq "W"){$h=$h*-1;} else {$h=$h+1;} push(@harray, $h);
			$v=substr $tile,5,2; $v=$v*1; $vsign=substr $tile,7,1; if ($vsign eq "S"){$v=$v*-1;} else {$v=$v+1;} push(@varray, $v);
			}
			@hsort = sort { $a <=> $b } @harray; $hmin = shift(@hsort); $hmax = pop(@hsort);
			@vsort = sort { $a <=> $b } @varray; $vmin = shift(@vsort); $vmax = pop(@vsort);
			$ulx=$hmin-1; $lry=$vmin-1; $lrx=$hmax; $uly=$vmax;
			$xsize=($lrx-$ulx)*$tileside; $ysize=($uly-$lry)*$tileside;
			$ulxcoord=$ulx; $ulycoord=$uly;
			#process tilelist
			open (OUT, ">aggregate_temp.txt");
			for $tile (@tilelist) {
			if (-e "$source/$name\_$year\_$tile.tif"){
			$h=substr $tile,0,3; $h=$h*1; $hsign=substr $tile,3,1; if ($hsign eq "W"){$h=$h*-1;} else {$h=$h+1;}
			$v=substr $tile,5,2; $v=$v*1; $vsign=substr $tile,7,1; if ($vsign eq "S"){$v=$v*-1;} else {$v=$v+1;}
			$x=($h-$hmin)*$tileside-2; $y=($vmax-$v)*$tileside-2;
			print OUT "$source/$name\_$year\_$tile.tif,$x,$y\n";
			}
			}
			close (OUT);
			readpipe "C:/GLAD_1.1/classification_aggregate.exe aggregate_temp.txt temp.img $xsize $ysize $ulxcoord $ulycoord >nul 2>&1";
			system "del aggregate_temp.txt";
			#build OVR
			open (BAT, ">temp_batch.bat");
			print BAT "gdal_translate -of GTiff -co COMPRESS=LZW -co BIGTIFF=IF_SAFER temp.img $outname.tif >nul 2>&1\n";
			print BAT "gdaladdo -ro $outname.tif 2 4 8 16 32 64 128 256 --config COMPRESS_OVERVIEW LZW >nul 2>&1\n";
			close (BAT);
			readpipe "\"$ogr\" temp_batch.bat >nul 2>&1";
			system "del temp.* >nul 2>&1";
			system "del temp_batch.bat >nul 2>&1";
			}
		}
	}
}
