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