next up previous contents
Next: About this document ... Up: SARAGO_REP Previous: Bibliography   Contents


APPENDIXES

#!/usr/bin/perl

# Routine di elaborazione dati SARAGO per creazione VOLUMI di ACQUA
# con salinita' determinate
#
# interfaccia un database MySQL e routines GMT
#
# Richiede una poligonale di inserimento dati
# e una superficie sintetica di fondo per effettuare i calcoli
# nello spazio indagato dal SARAGO
#

use Tplib;
use Getopt::Long;
use POSIX;

@optl = ("cond:s", "sal:s", "date:s", "range:s", "lines:s",
   "gen:s","var:s","grid","filter","pdf");
GetOptions @optl;

if ($opt_var eq "") {$opt_var="sal";}
#print "VAR: $opt_var\n";
$SAL_MIN = 10;
$SAL_MAX = 39.0;
$VAR{sal} = "5 30.0 38.0";
$VAR{temp} = "4 5 15";

if ($opt_sal)
{
   ($SAL_MIN,$SAL_MAX) = split ('/',$opt_sal);
   $VAR{sal} = "5 ${SAL_MIN} ${SAL_MAX}";
}

#$VAR{dens} = "5 15 30";
#$VAR{oxy} = "6 0 220";
#$VAR{cla} = "7 0 20";
#$VAR{prod} = "8 0 20";

$I="-I0.8m/1m";
$I="-I0.75m";
$I="-I0.25m/0.25m";
system "gmtset D_FORMAT %.10f";
system "rm -f ${opt_var}_vol.dat";

sub vmin {
 my (@D)=@_;
 my ($min) = 10e9;
 for ($d=0;$d<=$#D;$d++) {
  if ($D[$d]<$min) {$min=$D[$d];}
 }
 return $min;
}

sub vmax {
 my (@D)=@_;
 my ($max) = -10e9;
 for ($d=0;$d<=$#D;$d++) {
  if ($D[$d]>$max) {
     $max=$D[$d];
  }
 }
 return $max;
}
sub fill_volumes {
   my ($val)= @_;

}


PrjPar (1);
system "gmtset MEASURE_UNIT cm";

if ($opt_gen)
{
# data una linea generatrice AB e una distanza BC in m a sinistra (punto A a W)
# calcola le coordinate della linea parallela CD a distanza BC
# parallelepipedo che contiene i dati SARAGO e i campionamenti biologici
   print STDERR "Computing boxes ... " if $opt_v;
   $GEN = $opt_gen;
   $GEN =~ s/\://g;	
   ($lonA,$latA,$lonB,$latB,$dist) = split ('/',$GEN);
        $lonA = dmhd60($lonA);
        $lonB = dmhd60($lonB);
        $latA = dmhd60($latA);
        $latB = dmhd60($latB);
        ($D,$AZ12,$AZ21)= GeoPol($latA*$DTOR,$lonA*$DTOR,
   	$latB*$DTOR,$lonB*$DTOR);
        ($latC,$lonC,$az_back) = PolGeo ($latB*$DTOR,$lonB*$DTOR,
           			 $dist,$AZ12-90*$DTOR);
        ($latD,$lonD,$az_back) = PolGeo ($latA*$DTOR,$lonA*$DTOR,
           			 $dist,$AZ12-90*$DTOR);
   $latC *= $RTOD;
   $lonC *= $RTOD;
   $latD *= $RTOD;
   $lonD *= $RTOD;
# scrive su file per elaborazioni successive
   open POLY, ">POLY.DAT";
   print POLY "$lonA $latA\n$lonB $latB\n";
   print POLY "$lonC $latC\n$lonD $latD\n$lonA $latA\n";
        close POLY;
# calcola le coordinate del rettangolo che inscrive il parallelepipedo dati
# arrotondandole a 1' WSEN
   open BOX, ">BOX.DAT";
   $lonmin = vmin ($lonA,$lonB,$lonC,$lonD);
   $latmin = vmin ($latA,$latB,$latC,$latD);
   $lonmax = vmax ($lonA,$lonB,$lonC,$lonD);
   $latmax = vmax ($latA,$latB,$latC,$latD);
   $lonmin = sprintf "%.10f", dmhd60(int floor(d60dmh($lonmin)));
   $lonmax = sprintf "%.10f", dmhd60(int  ceil(d60dmh($lonmax)));
   $latmin = sprintf "%.10f", dmhd60(int floor(d60dmh($latmin)));
   $latmax = sprintf "%.10f", dmhd60(int  ceil(d60dmh($latmax)));
   print BOX ">\n";
   print BOX "$lonmin $latmin\n$lonmin $latmax\n";
   print BOX "$lonmax $latmax\n$lonmax $latmin\n$lonmin $latmin\n";
   close BOX;
   $R="-R${lonmin}/${lonmax}/${latmin}/${latmax}";
   print STDERR "Done with boxes\n" if $opt_v;
}

# Calcola la griglia di mask del box dati
if ($opt_grid)
{

}
($DATE_MIN,$DATE_MAX) = split ('/',$opt_date);

# recupera dati da database SARAGO MySQL
if ($opt_grid)
{
  print STDERR "Starting mask ... " if $opt_v;
  system "grdmask POLY.DAT -Gmask.grd -NNaN/1/1 $I $R";
  print STDERR "Done with mask \n" if $opt_v;
# Calcola la griglia di fondo sintetico DATI SARAGO sommando 0.5 m
  print STDERR "Starting synthetic bottom ... " if $opt_v;
  $command  = "awk '{print \$1,\$2,\$3+0.5}' sint_bat/SINT_BAT.DAT |";
  $command .= " surface -Ggrd/temp.grd $I $R -T0.05";
#  print "$command \n";
  system "$command";
# filtra ...
  system "grdfilter grd/temp.grd $R $I -D1 -Fb15 -Ggrd/filt.grd";
  system "grdmath grd/filt.grd mask.grd OR = grd/sint_bat.grd";
  print "Done with Synthetic bottom\n" if $opt_v;
  print STDERR "Start reading data ...\n" if ($opt_v);
  @lines = split (' ',$opt_lines);
  $i=0;
  for $l (@lines)
  {
   $QUERY  = "select LON,LAT,CAM_ID, LON,LAT,DEPTH,TEMP,SAL,";
   $QUERY .= "DATE_ACQ,TIME_ACQ,TZ,NFIX,LINE_ID ";
   $QUERY .= "from NAV where (LINE_ID = \'${l}\' ";
   $QUERY .="and DATE_ACQ>=\'${DATE_MIN}\' and DATE_ACQ<=\'${DATE_MAX}\'";
   $QUERY .= "and ( LAT >= ${latmin} and LAT<=${latmax}) ";
   $QUERY = $QUERY . "and (".$opt_cond.") " if $opt_cond;
   $QUERY .= "and ( LON >= ${lonmin} and LON<=${lonmax}) ";
   $QUERY .= "and SAL<=${SAL_MAX} and SAL>= ${SAL_MIN})";
   print STDERR "$QUERY\n" if ($opt_v);
   
   open SARAGO, "mysql_query localhost SARAGO \"${QUERY}\"|";
   while (<SARAGO>)
   {
   	chop;
   	($lat,$lon,@data) = split (',');
                if (($lat != $lat_old) or ($lon != $lon_old))
   	{
                  $DATA[$i] = "@data";
                  $i++;
   	}
   	$lat_old = $lat;
   	$lon_old = $lon;
   }
   print STDERR "Done with $l ...\n" if ($opt_v);
  }
  if ($opt_v)
  {
   for ($i=0;$i<=$#DATA;$i++)
   {
   	print STDERR "$DATA[$i]\n";
   }
  }
}

$d_step=0.5;
$d_step = 1;

open V,">VOLUMI.DAT" or die "Cannot open VOLUMI.DAT";
open TXT,">SARAGO.DAT";

# loop di calcolo per fascie di profondita' 1 m
$VTOT=0;
$AREA = ( 334.207 * 462.98);  # 0.25 nm !!!
$INF=32;
$SUP=37;
for ($v=0;$v<=($SUP-$INF+1);$v++) {$V[$v]=0;};
print "DATI SARAGO - ";
print "VOLUMI DI ACQUA (KM3) NELLE VARIE FRAZIONI DI SALINITA\n";
@cruise = split ('/',`pwd`);
print "CRUISE: ",$cruise[$#cruise-1],"\nAREA: ",uc $cruise[$#cruise];
printf "POLIGONO: (A) %.5f %.5f (B) %.5f %.5f (C) %.5f %.5f (D) %.5f %.5f\n",
   $lonA,$latA,$lonB,$latB,$lonC,$latC,$lonD,$latD,$lonA,$latA;
print "\nDEPTH \t  <${INF}\t";
for $i ($INF..$SUP-1) {
   print  "${i}-",$i+1,"\t";
} print "  >${SUP}  VTOT\n";

#
# per ogni fascia di profondita' fino a 35 m
#

for ($i=0;$i<35;$i+=$d_step)
 {
   $out=sprintf "%002d_%002d.dat",$i,$i+$d_step;
   $grd=sprintf "%002d_%002d.%s.grd",$i,$i+$d_step,$opt_var;
   ($N,$low,$high) = split (' ',$VAR{"sal"});
   if ($opt_grid) {
       open OUT,
         "|blockmean $R $I | surface -Gtemp.grd $R $I -Ll$low -Lu$high -T0.05";
       for $j (0..$#DATA)
       {
     @D=split(' ',$DATA[$j]);
     if ($D[3]>=$i and $D[3]<($i+$d_step+($i==0?0.5:0.0)))
     {
        print OUT "$D[1] $D[2] $D[$N]\n";
        print TXT "$D[1] $D[2] $D[3] $D[$N] $D[$#D]\n" if ($opt_v);
     }
   }
   close OUT;
     $opt_filter = 1;
          if ($opt_filter) {
           system "grdfilter temp.grd $R $I -D1 -Fb10 -Gfilt_grd";
   	system "grdmath filt_grd mask.grd OR = grd/${grd}";
   	print STDERR "Filtered ... " if $opt_v;
          }
     print "Done with $grd ...\n" if $opt_v;
   }

#  PLOT DATI SALINITA'


#  Clip con fondo sintetico DATI SARAGO


    open IN, "grd2xyz grd/$grd -V | grdtrack -Ggrd/sint_bat.grd | "
   	or die "Cannot open grdtrack ...";
    print "Start clipping $grd .." if $opt_v;
    $fname = "grd/".$grd;
    $fname =~ s/\.grd//;
    $fname .= "_masked.grd";
    $low = $i;
    $high = $i + $d_step;
    $depth = $high-0.5;
    open OUT , ">tempgrd";
    while (<IN>) {
   chop;
   ($lon,$lat,$val, $dsampled) = split;
   if ( ($depth>$dsampled) or ($dsampled =~/nan/i) ) {$val="NaN";}
    print OUT "$lon $lat $val\n";

   if ($val ne "NaN")
    {
           if ($val <= $SUP)
             {
               if  ($val <= $INF)
                {
                  $V[0]+=$AREA;
                } else {
               for $j ($INF..$SUP-1)
                     {
            		if ($val> $j  and $val <= $j+1)
            		 {
      	    		$V[$j+1-$INF]+=$AREA;
      	    		last;
            	         }
                }
             }
           } else {
               $V[($SUP-$INF)+1]+=$AREA;
           }
         }
     }
#
#   NITRATO = -2.72552 * SAL + 104
#   NITRITO = -0.242706 * SAL + 9.69066
#   SiO2 = (126.66 * SAL * SAL) - 7.53 * SAL + 0.112
#$NITRATO_INF34 = ((-2.72552 * 33.5 + 104) / 1.e6) * 62 ;

     close IN;
     close OUT;
     $VTOT = 0;
     printf "%002d-%002d ", $i,$i+1;
     for $j (0..$#V)
      {
    printf "%.5f ", $V[$j]/1.e9;
    $VTOT+=$V[$j]/1.e9;
      } printf "%.3f\n", $VTOT;
#
#     Costruisce griglia con maschera fondo sintetico
     system "xyz2grd tempgrd $R $I -G$fname" if $opt_grid;
     print STDERR "done with $fname .. \n" if ($opt_v and $opt_grid);
#
#     plot contours
     $opt_plot = 1;
     if ($opt_plot)
      {
         $J = "-JM0/44/15";
         $B = "-Bg1ma5mWSen";
         system "gmtset ANOT_FONT_SIZE 12";
         $B1= "-Bg1ma5mWSen";
         $ps = sprintf "sal-%002d-%002d.ps",$i,$i+1;
         $pdf = sprintf "sal-%002d-%002d.pdf",$i,$i+1;

         system "gmtset D_FORMAT %.2f";
         system "pscoast $R $J -W1 -Bg1ma5mWSen -Df -K -P   > ${ps}";
         system "grdcontour  ${fname} $R $J -C0.1 -A0.5 -K -O  >> ${ps} ";
#    system "psxy SARAGO.DAT -W1/255/0/0 -Sc0.1 $R $J -O -K >> ${ps} ";
         system "psxy POLY.DAT $R $J -M -O  -K >> ${ps}";
         open SC,">SC.DAT";
         print SC $lonmin+dmhd60(0.5)," ",$latmin+dmhd60(0.5),
               " 12 0 1 1 SAL ",$i,"-",$i+1," m ",$cruise[$#cruise-1]," ",
       	uc $cruise[$#cruise];
         close SC;
         system "pstext SC.DAT $R $J -K -O >> $ps ";
         system "echo \"showpage\" >> $ps ";
         system "convert $ps $pdf" if ($opt_pdf);
         system "gmtset D_FORMAT %.10f";
      }

#
#   dati volume  xyz binari
#      system "grd2xyz  $fname -Zf  >> ${opt_var}_vol.dat";

   
 }



G.Bortoluzzi 2001-07-16