#======================================================================
# L A D C P P R O C . L O A D C T D
# doc: Thu Dec 9 18:39:01 2010
# dlm: Sun Sep 14 21:20:12 2014
# (c) 2010 A.M. Thurnherr
# uE-Info: 37 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Dec 9, 2010: - exported from LADCPproc
# - added support for ASCII files
# Dec 16, 2010: - BUG cnv read did not work any more
# Jan 10, 2011: - added code to skip ANTS header
# Jan 22, 2011: - adapted to new -g
# Jul 15, 2011: - added $CTD{first_elapsed}
# Feb 5, 2012: - BUG: ASCII file did not deal with leading spaces correctly
# Apr 11, 2012: - BUG: ASCII file did not handle nan latlons correctly
# Apr 17, 2012: - fiddled
# Oct 19, 2012: - BUG: support for $CTD{first_elapsed} had not been implemented
# for binary files
# - BUG: CTD_badval had not been considered when setting $CTD{first_elapsed}
# - BUG: CNV format error was not detected correctly any more
# Jan 8, 2013: - added CTD_ASCII_header_lines
# Jun 25, 2013: - adapted to :: %PARAM convention
# Sep 25, 2013: - renamed "std" %PARAMs %lat, %lon, %ITS to conform to new convention
# - added support for carry-through of lat/lon info
# Nov 11, 2013: - BUG: lat/lon did not work any more?!?
sub readCTD_ASCII($$)
{
my($fn,$dtaR) = @_;
croak("$fn: unknown pressure field\n") unless defined($CTD_ASCII_press_field);
croak("$fn: unknown temperature field\n") unless defined($CTD_ASCII_temp_field);
croak("$fn: unknown salinity field\n") unless defined($CTD_ASCII_salin_field);
unless (numberp($dtaR->{stn_lat})) {
croak("$fn: unknown latitude field\n") unless defined($CTD_ASCII_lat_field);
croak("$fn: unknown longitude field\n") unless defined($CTD_ASCII_lon_field);
}
$CTD_ASCII_badval = 9e99 unless defined($CTD_ASCII_badval);
open(F,$fn) || croak("$fn: $!\n");
my($sumLat,$sumLon); my($nPos) = 0;
my($ds);
my($skip) = $CTD_ASCII_header_lines;
while (chomp($ds = <F>)) {
if (defined($CTD_ASCII_header_lines)) { # fixed header
next if ($skip-- > 0);
} else { # comments beginning with # allowed
next if ($ds =~ /^#/);
}
$ds =~ s/^\s+//; # strip leading spaces
my(@rec) = split('\s+',$ds);
push(@{$dtaR->{press}},($rec[$CTD_ASCII_press_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_press_field-1]);
push(@{$dtaR->{temp}}, ($rec[$CTD_ASCII_temp_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_temp_field-1]);
push(@{$dtaR->{salin}},($rec[$CTD_ASCII_salin_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_salin_field-1]);
if (defined($CTD_ASCII_lat_field) &&
numberp($rec[$CTD_ASCII_lat_field-1]) &&
$rec[$CTD_ASCII_lat_field-1] != $CTD_ASCII_badval) {
push(@{$dtaR->{lat}},$rec[$CTD_ASCII_lat_field-1]);
push(@{$dtaR->{lon}},$rec[$CTD_ASCII_lon_field-1]);
$nPos++;
$sumLat += $rec[$CTD_ASCII_lat_field-1];
$sumLon += $rec[$CTD_ASCII_lon_field-1];
} else {
push(@{$dtaR->{lat}},nan);
push(@{$dtaR->{lon}},nan);
}
}
close(F);
if ($nPos > 0) {
$dtaR->{stn_lat} = $sumLat / $nPos;
$dtaR->{stn_lon} = $sumLon / $nPos;
}
$dtaR->{sampint} = 1 / $CTD_ASCII_sampfreq;
}
sub readCTD_CNV($$)
{
my($fn,$dtaR) = @_;
my($CTD_nrecs,$CTD_nfields,$pressF,$tempF,$salinF,$elapsedF);
my($CTD_badval,$CTD_file_type);
open(F,$fn) || croak("$fn: $!\n");
while (1) { # parse header
my($hdr);
chomp($hdr = <F>);
croak(" unexpected EOF (format error)!\n") unless defined($hdr);
$hdr =~ s/\r*$//;
last if ($hdr eq '*END*');
$CTD_nfields = $',next if ($hdr =~ /nquan = /); # Layout
$CTD_nrecs = $',next if ($hdr =~ /nvalues = /);
$elapsedF = $1,next if ($hdr =~ /name (\d+) = timeS:/);
$pressF = $1,next if ($hdr =~ /name (\d+) = prDM:/);
if ($opt_2) {
$tempF = $1,next if ($hdr =~ /name (\d+) = t190C:/);
$salinF = $1,next if ($hdr =~ /name (\d+) = sal11:/);
} else {
$tempF = $1,next if ($hdr =~ /name (\d+) = t090C:/);
$salinF = $1,next if ($hdr =~ /name (\d+) = sal00:/);
}
$latF = $1,next if ($hdr =~ /name (\d+) = latitude:/);
$lonF = $1,next if ($hdr =~ /name (\d+) = longitude:/);
&antsAddParams('LADCPproc::CTD_start_time',$1),next # selected metadata
if ($hdr =~ /start_time = (.*)/);
&antsAddParams('LADCPproc::CTD_station',$1),next
if ($hdr =~ /Station\s*:\s*(.*)/);
&antsAddParams('LADCPproc::ship',$1),next
if ($hdr =~ /Ship\s*:\s*(.*)\s*$/);
&antsAddParams('LADCPproc::cruise',$1),next
if ($hdr =~ /Cruise\s*:\s*(.*)\s*$/);
&antsAddParams('LADCPproc::CTD_time',$1),next
if ($hdr =~ /Time\s*:\s*(.*)/);
&antsAddParams('LADCPproc::CTD_date',$1),next
if ($hdr =~ /Date\s*:\s*(.*)/);
if ($hdr =~ /Latitude\s*[=:]\s*/) {
($deg,$min,$NS) = split(/ /,$');
$dtaR->{stn_lat} = $deg + $min/60;
$dtaR->{stn_lat} *= -1 if ($NS eq 'S');
next;
}
if ($hdr =~ /Longitude\s*[=:]\s*/) {
($deg,$min,$EW) = split(/ /,$');
$dtaR->{stn_lon} = $deg + $min/60;
$dtaR->{stn_lon} *= -1 if ($EW eq 'W');
next;
}
if ($hdr =~ /interval = seconds: /) {
$dtaR->{sampint} = 1*$';
next;
}
$CTD_badval = $',next
if ($hdr =~ /bad_flag = /);
$CTD_file_type = $',next
if ($hdr =~ /file_type = /);
}
croak("$CTD_file: cannot determine CTD file layout\n")
unless ($CTD_nfields && $CTD_nrecs);
croak("$CTD_file: cannot determine missing value\n")
unless defined($CTD_badval);
croak("$CTD_file: not a CTD time series file\n")
unless ($dtaR->{sampint});
croak("$CTD_file: no pressure field\n")
unless defined($pressF);
croak("$CTD_file: no suitable temperature field\n")
unless defined($tempF);
croak("$CTD_file: no suitable salinity field\n")
unless defined($salinF);
if ($CTD_file_type eq 'ascii') {
while (1) {
last unless (my(@rec) = &antsFileIn(F));
$dtaR->{first_elapsed} = $rec[$elapsedF]
if !defined($dtaR->{first_elapsed}) && defined($elapsedF) && $rec[$elapsedF]!=$CTD_badval;
push(@{$dtaR->{press}},($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
push(@{$dtaR->{temp}}, ($rec[$tempF] == $CTD_badval) ? nan : $rec[$tempF]);
push(@{$dtaR->{salin}},($rec[$salinF] == $CTD_badval) ? nan : $rec[$salinF]);
push(@{$dtaR->{lat}},(!defined($latF) || ($rec[$latF] == $CTD_badval)) ? nan : $rec[$latF]);
push(@{$dtaR->{lon}},(!defined($lonF) || ($rec[$lonF] == $CTD_badval)) ? nan : $rec[$lonF]);
}
} elsif ($CTD_file_type eq 'binary') {
my($fbits) = 8 * length(pack('f',0));
croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
unless ($fbits == 32);
croak("$fn: can't read binary data\n")
unless (read(F,my($dta),4*$CTD_nfields*$CTD_nrecs) == 4*$CTD_nfields*$CTD_nrecs);
print(STDERR "$fn: WARNING: extraneous data at EOF\n") unless eof(F);
$dta = pack('V*',unpack('N*',$dta)) # big-endian CPU
if (unpack('h*', pack('s', 1)) =~ /01/); # c.f. perlport(1)
my(@dta) = unpack("f*",$dta);
for (my($r)=0; $r<$CTD_nrecs; $r++) {
$dtaR->{first_elapsed} = $dta[$r*$CTD_nfields+$elapsedF]
if !defined($dtaR->{first_elapsed}) && defined($elapsedF) && $dta[$r*$CTD_nfields+$elapsedF]!=$CTD_badval;
push(@{$dtaR->{press}},($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
push(@{$dtaR->{temp}}, ($dta[$r*$CTD_nfields+$tempF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$tempF]);
push(@{$dtaR->{salin}},($dta[$r*$CTD_nfields+$salinF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$salinF]);
push(@{$dtaR->{lat}},(!defined($latF) || ($dta[$r*$CTD_nfields+$latF] == $CTD_badval)) ? nan : $dta[$r*$CTD_nfields+$latF]);
push(@{$dtaR->{lon}},(!defined($lonF) || ($dta[$r*$CTD_nfields+$lonF] == $CTD_badval)) ? nan : $dta[$r*$CTD_nfields+$lonF]);
}
} else {
croak("$fn: unknown CTD file type $CTD_file_type\n");
}
close(F);
}
sub readCTD($$)
{
my($fn,$dtaR) = @_;
if (defined($CTD_ASCII_sampfreq)) {
readCTD_ASCII($fn,$dtaR);
} else {
readCTD_CNV($fn,$dtaR);
}
croak("$0: unknown latitude\n") unless defined($dtaR->{stn_lat});
&antsAddParams('LADCPproc::CTD_lat',$dtaR->{stn_lat});
croak("$0: unknown longitude\n") unless defined($dtaR->{stn_lon});
&antsAddParams('LADCPproc::CTD_lon',$dtaR->{stn_lon});
&antsAddParams('LADCPproc::CTD_sampfreq',1/$dtaR->{sampint});
&antsAddParams('LADCPproc::CTD_ITS',$P{ITS} = 90);
}
1;