LADCPproc.loadCTD
author A.M. Thurnherr <athurnherr@yahoo.com>
Sun, 27 Jul 2014 16:28:57 -0400
changeset 30 8697ba5a88ec
parent 24 40756f8aff85
child 29 f72cd642972c
permissions -rw-r--r--
before beginning to add code to calculate package velocity

#======================================================================
#                    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: Thu Nov 21 10:14:37 2013
#                    (c) 2010 A.M. Thurnherr
#                    uE-Info: 170 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;