2011_07_10/LADCPproc.loadCTD
changeset 10 196a179304ee
parent 9 48b2d27aaebf
child 11 d0af7f7aa23b
equal deleted inserted replaced
9:48b2d27aaebf 10:196a179304ee
     1 #======================================================================
       
     2 #                    L A D C P P R O C . L O A D C T D 
       
     3 #                    doc: Thu Dec  9 18:39:01 2010
       
     4 #                    dlm: Sat Jan 22 21:40:12 2011
       
     5 #                    (c) 2010 A.M. Thurnherr
       
     6 #                    uE-Info: 174 0 NIL 0 0 72 2 2 4 NIL ofnI
       
     7 #======================================================================
       
     8 
       
     9 # HISTORY:
       
    10 #	Dec  9, 2010: - exported from LADCPproc
       
    11 #				  - added support for ASCII files
       
    12 #	Dec 16, 2010: - BUG cnv read did not work any more
       
    13 #	Jan 10, 2011: - added code to skip ANTS header
       
    14 #	Jan 22, 2011: - adapted to new -g
       
    15 
       
    16 sub readCTD_ASCII($$)
       
    17 {
       
    18 	my($fn,$dtaR) = @_;
       
    19 
       
    20 	croak("$fn: unknown pressure field\n")    unless defined($CTD_ASCII_press_field);
       
    21 	croak("$fn: unknown temperature field\n") unless defined($CTD_ASCII_temp_field);
       
    22 	croak("$fn: unknown salinity field\n")    unless defined($CTD_ASCII_salin_field);
       
    23 	unless (numberp($dtaR->{lat})) {
       
    24 		croak("$fn: unknown latitude field\n")    unless defined($CTD_ASCII_lat_field);
       
    25 		croak("$fn: unknown longitude field\n")   unless defined($CTD_ASCII_lon_field);
       
    26 	}
       
    27 
       
    28 	$CTD_ASCII_badval = 9e99 unless defined($CTD_ASCII_badval);
       
    29 
       
    30 	open(F,$fn) || croak("$fn: $!\n");
       
    31 	my($sumLat,$sumLon); my($nPos) = 0;
       
    32 	my($ds);
       
    33 	while (chomp($ds = <F>)) {
       
    34 		next if ($ds =~ /^#/);
       
    35 		my(@rec) = split('\s+',$ds);
       
    36 		push(@{$dtaR->{press}},($rec[$CTD_ASCII_press_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_press_field-1]);
       
    37 		push(@{$dtaR->{temp}}, ($rec[$CTD_ASCII_temp_field-1]  == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_temp_field-1]);
       
    38 		push(@{$dtaR->{salin}},($rec[$CTD_ASCII_salin_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_salin_field-1]);
       
    39 		unless (!defined($CTD_ASCII_lat_field) || $rec[$CTD_ASCII_lat_field-1] == $CTD_ASCII_badval) {
       
    40 			$nPos++;
       
    41 			$sumLat += $rec[$CTD_ASCII_lat_field-1];
       
    42 			$sumLon += $rec[$CTD_ASCII_lon_field-1];
       
    43 		}
       
    44 	}
       
    45 	close(F);
       
    46 	
       
    47 	if ($nPos > 0) {
       
    48 		$dtaR->{lat} = $sumLat / $nPos;
       
    49 		$dtaR->{lon} = $sumLon / $nPos;
       
    50 	}
       
    51 
       
    52 	$dtaR->{sampint} = 1 / $CTD_ASCII_sampfreq;
       
    53 }
       
    54 
       
    55 sub readCTD_CNV($$)
       
    56 {
       
    57 	my($fn,$dtaR) = @_;
       
    58 	my($CTD_nrecs,$CTD_nfields,$pressF,$tempF,$salinF);
       
    59 	my($CTD_badval,$CTD_file_type);
       
    60 
       
    61 	open(F,$fn) || croak("$fn: $!\n");
       
    62 	while (1) { 														# parse header
       
    63 		my($hdr);
       
    64 		chomp($hdr = <F>);
       
    65 		$hdr =~ s/\r*$//;
       
    66 		croak("$0: unexpected EOF (format error)\n") unless defined($hdr);
       
    67 		last if ($hdr eq '*END*');
       
    68 	    
       
    69 		$CTD_nfields = $',next if ($hdr =~ /nquan = /); 				# Layout
       
    70 		$CTD_nrecs = $',next if ($hdr =~ /nvalues = /);
       
    71 		$pressF = $1,next if ($hdr =~ /name (\d+) = prDM:/);
       
    72 		if ($opt_2) {
       
    73 			$tempF	= $1,next if ($hdr =~ /name (\d+) = t190C:/);
       
    74 			$salinF = $1,next if ($hdr =~ /name (\d+) = sal11:/);
       
    75 		} else {
       
    76 			$tempF	= $1,next if ($hdr =~ /name (\d+) = t090C:/);
       
    77 			$salinF = $1,next if ($hdr =~ /name (\d+) = sal00:/);
       
    78 		}
       
    79 	
       
    80 		&antsAddParams('start_time',$1),next							# selected metadata
       
    81 			if ($hdr =~ /start_time = (.*)/);
       
    82 	
       
    83 		&antsAddParams('station',$1),next
       
    84 			if ($hdr =~ /Station\s*:\s*(.*)/);
       
    85 		&antsAddParams('ship',$1),next
       
    86 			if ($hdr =~ /Ship\s*:\s*(.*)\s*$/);
       
    87 		&antsAddParams('cruise',$1),next
       
    88 			if ($hdr =~ /Cruise\s*:\s*(.*)\s*$/);
       
    89 		&antsAddParams('time',$1),next
       
    90 			if ($hdr =~ /Time\s*:\s*(.*)/);
       
    91 		&antsAddParams('date',$1),next
       
    92 			if ($hdr =~ /Date\s*:\s*(.*)/);
       
    93 	
       
    94 		if ($hdr =~ /Latitude\s*[=:]\s*/) {
       
    95 			($deg,$min,$NS) = split(/ /,$');
       
    96 			$dtaR->{lat} = $deg + $min/60;
       
    97 			$dtaR->{lat} *= -1 if ($NS eq 'S');
       
    98 			next;
       
    99 		}
       
   100 		if ($hdr =~ /Longitude\s*[=:]\s*/) {
       
   101 			($deg,$min,$EW) = split(/ /,$');
       
   102 			$dtaR->{lon} = $deg + $min/60;
       
   103 			$dtaR->{lon} *= -1 if ($EW eq 'W');
       
   104 			next;
       
   105 		}
       
   106 	    
       
   107 		if ($hdr =~ /interval = seconds: /) {
       
   108 			$dtaR->{sampint} = 1*$';
       
   109 			next;
       
   110 		}
       
   111 	    
       
   112 		$CTD_badval = $',next
       
   113 			if ($hdr =~ /bad_flag = /); 
       
   114 		$CTD_file_type = $',next
       
   115 			if ($hdr =~ /file_type = /);    
       
   116 	}
       
   117 	
       
   118 	croak("$CTD_file: cannot determine CTD file layout\n")
       
   119 		unless ($CTD_nfields && $CTD_nrecs);
       
   120 	croak("$CTD_file: cannot determine missing value\n")
       
   121 		unless defined($CTD_badval);
       
   122 	croak("$CTD_file: not a CTD time series file\n")
       
   123 		unless ($dtaR->{sampint});
       
   124 	croak("$CTD_file: no pressure field\n")
       
   125 		unless defined($pressF);
       
   126 	croak("$CTD_file: no suitable temperature field\n")
       
   127 		unless defined($tempF);
       
   128 	croak("$CTD_file: no suitable salinity field\n")
       
   129 		unless defined($salinF);
       
   130 	
       
   131 	if ($CTD_file_type eq 'ascii') {
       
   132 		while (1) {
       
   133 			last unless (my(@rec) = &antsFileIn(F));
       
   134 			push(@{$dtaR->{press}},($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
       
   135 			push(@{$dtaR->{temp}}, ($rec[$tempF]  == $CTD_badval) ? nan : $rec[$tempF]);
       
   136 			push(@{$dtaR->{salin}},($rec[$salinF] == $CTD_badval) ? nan : $rec[$salinF]);
       
   137 		}
       
   138 	} elsif ($CTD_file_type eq 'binary') {
       
   139 	
       
   140 		my($fbits) = 8 * length(pack('f',0));
       
   141 		croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
       
   142 			unless ($fbits == 32);  
       
   143 		    
       
   144 		croak("$fn: can't read binary data\n")
       
   145 			unless (read(F,my($dta),4*$CTD_nfields*$CTD_nrecs) == 4*$CTD_nfields*$CTD_nrecs);
       
   146 		print(STDERR "$fn: WARNING: extraneous data at EOF\n") unless eof(F);
       
   147 	
       
   148 		$dta = pack('V*',unpack('N*',$dta)) 			# big-endian CPU
       
   149 			if (unpack('h*', pack('s', 1)) =~ /01/);	# c.f. perlport(1)
       
   150 	    
       
   151 		my(@dta) = unpack("f*",$dta);
       
   152 	
       
   153 		for (my($r)=0; $r<$CTD_nrecs; $r++) {
       
   154 			push(@{$dtaR->{press}},($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
       
   155 			push(@{$dtaR->{temp}}, ($dta[$r*$CTD_nfields+$tempF]  == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$tempF]);
       
   156 			push(@{$dtaR->{salin}},($dta[$r*$CTD_nfields+$salinF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$salinF]);
       
   157 		}
       
   158 	} else {
       
   159 		croak("$fn: unknown CTD file type $CTD_file_type\n");
       
   160 	}
       
   161 	close(F);
       
   162 }
       
   163 
       
   164 sub readCTD($$)
       
   165 {
       
   166 	my($fn,$dtaR) = @_;
       
   167 
       
   168 	if (defined($CTD_ASCII_sampfreq)) {
       
   169 		readCTD_ASCII($fn,$dtaR);
       
   170 	} else {
       
   171 		readCTD_CNV($fn,$dtaR);
       
   172 	}
       
   173 
       
   174 	croak("$0: unknown latitude\n") unless defined($dtaR->{lat});
       
   175 	&antsAddParams('lat',$dtaR->{lat});
       
   176 	croak("$0: unknown longitude\n") unless defined($dtaR->{lon});
       
   177 	&antsAddParams('lon',$dtaR->{lon});
       
   178 
       
   179 	&antsAddParams('CTD_sampfreq',1/$dtaR->{sampint});
       
   180 	&antsAddParams('ITS',$P{ITS} = 90);
       
   181 }
       
   182 
       
   183 1;