#!/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";
}