#!/usr/bin/perl
#compareDEMs.pl -- Compares two sorted xyz DEM files and reports the mean  
#topography deviation

use Getopt::Std;
getopts('a:b:');

$usage = "compareDEMs.pl -a <first DEM xyz file> -b <second DEM xyz file>";

if (! ($opt_a && $opt_b )) { die "\nUsage: $usage\n\n"; }

open(DEM1,"$opt_a") || die "Could not open track file.\n";

open(DEM2,"$opt_b") || die "Could not open data file.\n";

$counter = 0;
$elevation_sum =0;
$invalid_elevation = 0;

while(<DEM1>) {

  $line = $_;
  if ( $line =~ /^(\S+)\s+(\S+)\s+(\S+)/ ) {
   #Grab the first two columns of data

     $longitude = $1;
     $latitude = $2;
     $elevation = $3;

     $match = 0;
     while ( $match == 0 )  {
        $line2 = <DEM2>;
        if ( $line2 =~ /^(\S+)\s+(\S+)\s+(\S+)\s/ ) {
           $longitude2 = $1;
           $latitude2 = $2;
           $elevation2 = $3;


           #Check if the coordinates match

           if ( ( $longitude == $longitude2) && ( $latitude == $latitude2 ) ) {
              $match = 1;
              if (($elevation != $invalid_elevation) && ($elevation2 != $invalid_elevation )) {
                @elevation_difference[$counter] = abs($elevation - $elevation2);
                $elevation_sum = $elevation_sum + @elevation_difference[$counter]; 
                $counter++;
              }  

           }
        }
     }
  }
}

$temp_sum = 0;
$mean = $elevation_sum / $counter;

for ($x=0; $x < $counter; $x++) {
  $temp_sum = $temp_sum + (@elevation_difference[$x] - $mean)**2;
}
$variance = $temp_sum / $counter;
$standard_deviation = sqrt($variance);  
print "The mean topography deviation is $mean\n";
print "The standard deviation is $standard_deviation\n";