#!/usr/bin/perl

use POSIX qw(ceil floor);
#use IO::File;

##### Usage info/check and assignment of command line arguments

sub Usage {
  print STDERR <<END;

  Usage: mosaic_srtm.pl bounds dem_name [dir_destination] [DIR_data]
    Function:  Builds DEM for a specified area using SRTM data.
    Bounds must be in -Rw/e/s/n format.
    dir is the destination directory (default is current)

  *** Yuri Fialko, Sep 21, 2002 *** 
END
  exit 1;
}
@ARGV > 1 or Usage();

$bounds=shift;
$outdem=shift;
$dird=shift or chomp($dird=`pwd`);
$dir0=shift or $dir0 = "/winsar/DEM/SRTM/United_States/1arcsec/";
#$dir0=shift or $dir0 = "/net/winsar/DEM/SRTM/North_America/3arcsec/";
#$dir0 = "/winsar/DEM/SRTM/World/";
##$dird=$dird."/";
eval(`cd $dird`);
print STDERR "$dird \n";

#gets limits from $bounds
($bdash,$bR,$bcoords) = split (//,$bounds,3);
    ($bw,$be,$bs,$bn) = split /\//,$bcoords;

$lon1=floor($bw); 
$lon2=ceil($be)-1; 
$lat1=floor($bs); 
$lat2=ceil($bn)-1; 
# May  not work for stitching actoss the equator!

for $lat($lat1..$lat2){ 
 if ($lat <= 0){$LT="S"; $latc=-$lat}
 else{$LT="N"; $latc=$lat}
 $latc=$LT."$latc";
 if (!-e "$dird/lon.$latc.dem"){
  for $lon($lon1+1..$lon2){ 
   if ($lon == $lon1+1){
    $lon_noodle=$lon-1;
    if ($lon_noodle <= 0){$LN="W"; $lon_noodle=-$lon_noodle}
    else{$LN="E"}
    if ($lon_noodle < 10){$lon_noodle="00".$lon_noodle;}
    else { if ($lon_noodle < 100){$lon_noodle="0".$lon_noodle;}}
    $fname=$latc.$LN."$lon_noodle";
    $add_noodle=$dir0.$fname;
   }
   else{$add_noodle=$dird."/lon.$latc"}
   if ($lon <= 0){$LN="W"; $lonc=-$lon}
    else{$LN="E"; $lonc=$lon}
    if ($lonc < 10){$lonc="00".$lonc;}
    else { if ($lonc < 100){$lonc="0".$lonc;}}
    $lonc=$LN."$lonc";
    $fname=$latc.$lonc;
    $ffname=$dir0.$fname;
    $target = $dird."/lon.".$latc;
    print STDERR "Adding stamp $latc $lonc \n";
      print STDERR "SRTMpaste $ffname $add_noodle $target \n";
    `SRTMpaste $ffname $add_noodle $target`;
  }
 }
 if ($lon1 == $lon2){
  if ($lon1 <= 0){$LN="W"; $lonc=-$lon1}
   else{$LN="E"; $lonc=$lon1}
   $lonc=$LN."$lonc";
   $fname=$latc.$lonc;
   $ffname=$dir0.$fname;
print STDERR "cp $ffname* $dird \n"; 
    `cp $ffname* $dird`;
    `chmod 644 $fname*`;
    `mmv.pl $fname lon.$latc $dird`;
 }
 if ($lat1 == $lat2){
   $latp=$lat-1;
   $latp=$lastlat.abs($latp);
  `mmv.pl lon.$latc $outdem $dird`;
 }
 if ($lat > $lat1) {
   if ($lat == $lat1+1){
    $latp=$lat-1;
    $latp="lon.".$lastlat.abs($latp);
    $add_noodle=$latp}
   else{$add_noodle=$outdem}
   print STDERR "Stitching along $latc \n";
   print STDERR "SRTMpaste $dird/lon.$latc $dird/$add_noodle $dird/$outdem \n";
  `SRTMpaste $dird/lon.$latc $dird/$add_noodle $dird/$outdem`;
 }
   $lastlat=$LT;
}
open ORBIN, ">>$dird/$outdem.dem.rsc" or die "Can't write to $outdem.dem.rsc\n";
 print ORBIN "Z_OFFSET                0  \n";
 print ORBIN "Z_SCALE                 1  \n";
 print ORBIN "PROJECTION              LL \n";
 print ORBIN "RLOOKS                  1  \n";
 print ORBIN "ALOOKS                  1  \n";
close(ORBIN);
open ORBIN, ">>$dird/$outdem.dem.rsc.hst" or die "Can't write to $outdem.dem.rsc\n";
 print ORBIN "mosaic_srtm.pl $bounds $outdem  \n";
close(ORBIN);
`rm -f $dird/lon.*`;
print STDERR "Finished building DEM $outdem \n";