Revision Log

tower_parse_dev.perl

Revision 7 (checked in by cbc, 18 years ago)	--
Property svn:executable set to

Line
1	#!/afs/isis/pkg/perl-582/bin/perl
2	#tower_strip.perl
3	#Time-stamp: <2006-06-12 11:52:14 lstearns>
4	#
5	#Abstract:
6	# This script is designed to parse existing tower data for a given
7	# time range and save the output as a netcdf file containing the
8	# converted data values as well as ancillary information.
9	#
10	# ./tower_parse.perl -help
11	# For more information on input arguments
12	#
13	#Created: 2003.09?.?? maybe
14	#Author: Luke Stearns
15	#Modifications: started back up 2004.03.?? still going....
16	#
17	# Luke Stearns 2006.03.15 - going back to improve some of the comments
18	#
19
20
21	######################################################################
22	# Set up libraries and paths.
23
24	use lib '/opt/local/seacoos/bin/';
25	use lib '/opt/local/seacoos/bin/meta_data/';
26	use meta_data;
27	$bin_dir = '/opt/local/seacoos/bin/meta_data/';
28
29	#meta_data.pm is a module I wrote to access some data tables in
30	#/opt/local/seacoos/bin/meta_data/
31	#I intended these tables as a first step towards designing a relational
32	#database to store instrument configuration and deployment information.
33
34	#Development paths and libraries:
35	# use lib '/afs/isis.unc.edu/depts/marine/workspace/hseim/sablam4/'.
36	# 'lstearns/perl/test_nemo/';
37	# use lib '/afs/isis.unc.edu/depts/marine/workspace/hseim/sablam4/'.
38	# 'lstearns/perl/test_nemo/meta_data/';
39	# use meta_data;
40	# $bin_dir = '/afs/isis.unc.edu/depts/marine/workspace/hseim/sablam4/'.
41	# 'lstearns/perl/test_nemo/meta_data/';
42
43	#Other included modules:
44	use NetCDF;
45	use Date::Manip;
46	use Statistics::Basic::StdDev;
47	use Math::BigFloat;
48
49	######################################################################
50	# Keep track of the time it takes to run the script
51	$starttime = time;
52
53	######################################################################
54	# Handle input arguments, namely a date string and a path, if given
55	# I'd like to handle a variety of input arguments
56	&handle_input;
57
58	######################################################################
59	# Read information about variables from the tables: parsing info,
60	# conversions, and meta-data
61
62	#This is now handled by the 'use meta_data;' line above.
63
64	######################################################################
65	# Read in the data files and parse each line of data!
66	# This subroutine drives the bulk of the script!
67	if ($debug) {print "\nReading data files:\n";}
68	$first_one = 1;
69	&read_data;
70
71	######################################################################
72	# Finally, write the data to the NetCDF file. This is where it
73	# crashes more often than not.
74	#&test_data;
75	&write_netcdf;
76
77	######################################################################
78	# Let us know how long it took!
79	$cummtime = time - $starttime;
80	if ($debug) {
81	print "Closing NetCDF files:\nALL DONE!";
82	print "\nTotal time $cummtime sec\n";
83	}
84
85	if ($nemo) {
86	print "Attempting to scp file to nemo... (will require passwd)\n";
87	system("scp $out_data_dir$out_file_name nemo.isis.unc.edu:/seacoos/".
88	"data/nc-coos/roof/");
89	}
90
91	######################################################################
92	####################### BEGIN SUBROUTINES ############################
93	######################################################################
94
95	######################################################################
96	### Read Data ########################################################
97	sub read_data {
98
99	#Set the starting time index (will be incremented to 0 first)
100	$time_index = -1;
101	$file_index = -1;
102
103	#Set the fill value in case of missing data.
104	$fill_value = -9999;
105
106	#Get all of the filenames (including path) to the files in the
107	#directory of data that satisfy the matching pattern.
108	@tower_files = ();
109
110	#For each of the dates determined by the user input (&handle_input
111	#takes us to &handle_date that sets the array @dates) find the files
112	#for that date and store them in the array @tower_files
113	foreach $date_string (@dates) {
114	$data_sub_dir = $month_dirs{$date_string};
115	@t_files = glob("${d_data_dir}/$data_sub_dir/${tower_id}_${date_string}.tx");
116	push(@tower_files,@t_files);
117
118	}
119
120	if ($debug) {
121	print "... Found ".eval($#tower_files +1)." files.\n";
122	}
123
124	#I'm not sure where else to put this... I want to hang on to the
125	#sample loop time, and this is where I put it...
126	$expected_vars{time} = "SAMPLE_LOOP";
127
128	#Loop through the filenames loading, then reading each file...
129	foreach $file (@tower_files) {
130	$file_time = $file;
131	#$file_time =~m/([0-9_]+).txt/;
132	$file_time =~m/([0-9]{4})([0-9]{2})([0-9]{2})_([0-9]{2})([0-9]{2})[0-9]+.txt/;
133	$file_time_str = "$1/$2/$3 $4:$5:00";
134
135	#Convert to number of seconds since epoch
136	$file_u_date = &UnixDate(&ParseDate("$file_time_str UTC"),"%s");
137
138	#1970/01/01 00:00:00 as a datenum
139	$epoch_dn = 719529;
140
141	#Convert the seconds since epoch to a datenum
142	$file_datenum = $file_u_date/86400 + $epoch_dn;
143
144	$file_index++;
145	${$Vars{'file_time'}}[$file_index] = $file_datenum;
146
147	open(DATA_FILE,$file);
148	#Loop through all lines of data
149	while ($data_line=<DATA_FILE>) {
150
151	#Here's the work! Seperate all of the data and store it in
152	#a way that we can easily stuff it into the NetCDF files at
153	#the end of it all
154	&seperate_sensors;
155
156	} #Done reading this file
157	close(DATA_FILE);
158
159	} #Done looping through the files
160
161	#Just a thought here... perhaps all parsing and conversions should
162	#occur here, operating on the arrays of data, rather than for each
163	#time... That might also be a better way to handle derived
164	#variables like pir, density, etc. Maybe later.
165	}
166
167	######################################################################
168	### Test Data ########################################################
169	sub test_data {
170	print "Trying all keys of data array\n";
171	foreach $param (@vars) {
172	print "...$param,";
173	@temp_array = @{$Vars{$param}};
174	print "Array: $#temp_array, $temp_array[0]\n";
175	}
176
177	}
178	######################################################################
179	### Write Netcdf #####################################################
180	sub write_netcdf {
181
182	#I think I'd like to pull this out into another script. The point
183	#of this subroutine is to initialize and populate the NetCDF file
184	#that will contain the data.
185
186	####################################################################
187	# Create a NetCDF file to work with!
188	if ($debug) {
189	print "\nCreating NetCDF files:\n";
190	}
191	$nc_name = "${out_data_dir}$out_file_name";
192	$nc_sname = "${out_data_dir}$out_file_sname";
193
194	unless ($#platform_time>0) {
195	print "Exiting Now! \nFATAL ERROR:\n";
196	die " There must be at least one valid record!\n".
197	" Please try another date range!\n";
198	}
199
200
201	$ncid = NetCDF::create($nc_name,NetCDF::CLOBBER);
202	die "Couldn't open $nc_name\n" if $ncid<0;
203
204	$ncid_system = NetCDF::create($nc_sname,NetCDF::CLOBBER);
205	die "Couldn't open $nc_sname\n" if $ncid_system<0;
206
207
208	####################################################################
209	# Put Global Attributes...
210
211	if ($debug) {
212	print "... Putting global variables\n";
213	}
214
215	#I had this, but it puts them in random order, so I'll use the order
216	#direct from globals.table
217	foreach $global_att (@globals) {
218	$g_a_value = eval("$global_value{$global_att}");
219	$g_a_format = $global_format{$global_att};
220
221	unless ($g_a_value =~/N\/A/) {
222	if ($debug>2) {
223	print "$global_att is $g_a_value\n"
224	}
225	$attid = NetCDF::attput($ncid,NetCDF::GLOBAL,"$global_att",
226	eval('NetCDF::'.$g_a_format),
227	"$g_a_value");
228	$attid = NetCDF::attput($ncid_system,NetCDF::GLOBAL,"$global_att",
229	eval('NetCDF::'.$g_a_format),
230	"$g_a_value");
231	$attid != -1 \|\| die "\n Couldn't write global attribute $global_att\n";
232	}
233	}
234
235
236
237	####################################################################
238	# Define dimensions, fixed ones first
239	@{$Vars{LONGITUDE}} = ($tower_lon{$tower_id});
240	@{$Vars{LATITUDE}} = ($tower_lat{$tower_id});
241
242	#I did it this way because I was starting from Brians script (which
243	#was generally extremely helpful)
244	%dims = (
245	LONGITUDE => 1,
246	LATITUDE => 1,
247	TOWER_FILES => $#tower_files+1,
248	CHARACTERS => 50,
249	time => $#platform_time+1
250	);
251	%dim_names = (
252	'a' => 'LATITUDE',
253	'o' => 'LONGITUDE',
254	'f' => 'TOWER_FILES',
255	'c' => 'CHARACTERS',
256	't' => 'time'
257	);
258
259
260	foreach $dim (keys %dims) {
261	if ($dimensions{$dim}=~/[fc]/){
262	$ncid{$dim} = $ncid_system;
263	}
264	elsif ($dimensions{$dim}=~/[tao]/) {
265	$ncid{$dim} = $ncid;
266	}
267
268	unless ($dim =~/time/) {
269	$dimid{$dim} = NetCDF::dimdef($ncid{$dim},$dim,$dims{$dim});
270	die "Couldn't define $dim dimension\n" if $dimid{$dim}<0;
271	}
272	}
273
274	####################################################################
275	# Define time dimension
276	$dimid{time} = NetCDF::dimdef($ncid,"time",NetCDF::UNLIMITED);
277	die "Couldn't define record dimension\n" if $dimid{time}<0;
278
279
280	####################################################################
281	# Define variables and put attributes, fixed ones first
282	if ($debug) {
283	print "... Defining variables and putting attributes\n";
284	}
285	####################################################################
286	# Define variables and put attributes, now time dependant ones
287
288	#Dimension variables like lat and lon.
289	@all_attributes = ('long_name', 'standard_name', 'units',
290	'description', 'range', 'precision'
291	);
292
293	@zeros = (0,0,0); #Variables are dependant at most on 4 dimensions
294
295	#foreach $param (keys %Vars) {
296	foreach $param (@vars) {
297	if ($param =~/CLOCK/ \|\| $param =~/BBAD/) {
298	#We'll skip these for now...
299	}
300	elsif (defined $Vars{$param}) {
301	&vardef_attput;
302	}
303	}
304
305	####################################################################
306	# End Definition period!
307	$status = NetCDF::endef($ncid);
308	die "Couldn't end definition\n" if $status<0;
309	$status = NetCDF::endef($ncid_system);
310	die "Couldn't end definition\n" if $status<0;
311
312
313	####################################################################
314	# Start putting variables, first the fixed variables
315	if ($debug) {
316	print "... Putting variables:\n\n";
317	}
318	#foreach $param (keys %Vars) {
319	foreach $param (@vars) {
320	if ($param =~/CLOCK/ \|\| $param =~/BBAD/) {
321	#We'll skip these for now...
322	}
323	elsif (defined $Vars{$param}) {
324	&varput;
325	}
326	}
327
328	####################################################################
329	# All done, let's close up!
330	$status = NetCDF::close($ncid);
331	$status = NetCDF::close($ncid_system);
332
333	}
334
335	######################################################################
336	### Seperate Sensors #################################################
337	sub seperate_sensors {
338
339	#This subroutine parses each line of data in a data file: determines
340	#what the format of the line is, and sends the results onward to
341	#compute derived variables, parse complicated data strings, and
342	#convert to scientific units. The organization of this one is
343	#pretty confusing...
344
345	@data_line = split(', ',$data_line);
346	@va = ();
347	#If there is at least one comma in the data line it is of the usual
348	#data format (not ADCP!) and we'll treat it normally.
349
350	if ($#data_line) {
351	#Everything other than the CTD and ADCP should fit this...
352	if ($#data_line == 5) {
353	($time,$param_id,$type,$value,$units,$error_string) = @data_line;
354	}
355
356	#The CTD is still a special case becuase data is comma delimited.
357	elsif ($#data_line == 8) {
358	($time,$param_id,$type,$va[0],$va[1],$va[2],$va[3],$units,
359	$error_string) = @data_line;
360	$value = join(',',@va);
361	$value =~s/\s+//g;
362	}
363
364	#The CTD is still a special case becuase data is comma delimited.
365	elsif ($#data_line == 9) {
366	#print "9\n@data_line\n";
367
368	($time,$param_id,$type,$va[0],$va[1],$va[2],$va[3],$va[4],
369	$units,$error_string) = @data_line;
370	$value = join(',',@va);
371	$value =~s/\s+//g;
372	}
373
374	#The CTD is still a special case becuase data is comma delimited.
375	elsif ($#data_line == 10) {
376	#print "10\n@data_line\n";
377
378	($time,$param_id,$type,$va[0],$va[1],$va[2],$va[3],$va[4],
379	$va[5],$units,$error_string) = @data_line;
380	$value = join(',',@va);
381	$value =~s/\s+//g;
382	}
383
384	else {
385	#print "WARNING: Data line not parsed correctly!\n@data_line\n";
386	#This is so a failed parse doesn't end up going to the previous
387	#param_id
388	($time,$param_id,$type,@junk) = @data_line;
389	#print "!!!!$param_id\n";
390
391	#$param_id = 'failed';
392	$value = "$fill_value,$fill_value,$fill_value,$fill_value";
393
394	}
395	if($param_id=~m/CTD1/) {
396	#$press = 0;
397	#print "HEY! @data_line";
398	#print "$sal = $temp = $press = $cond = $fill_value\n";
399	}
400
401	#If it has failed the previous two options, than something is wrong.
402	#In fact, the second choice shouldn't ever occur, because Sara has
403	#implemented a change that puts quotes around data values that need
404	#special parsing.
405
406	#elsif ($debug>2) {
407	#print "Error suspected for data line:\n$data_line\n";
408	#}
409	$big_type = $type;
410	$big_type =~s/$[0-9]$//g;
411
412	#This is a first attempt at dealing with variables that contain
413	#a number of sub variables. For example the CTD output string
414	#contains temparature conductivity salinity and density.
415	#our (@sub_vars,%sub_vars);
416
417	@sub_vars = @{$sub_vars{$param_id}};
418
419	&compute_sub_vars;
420
421
422	#If we are at the start of a file, keep track of the settings
423	#listed, especially the param_id's that we expect to be in
424	#each record, so that we can insert fill values where needed.
425
426	if ($type=~/SETTING/ && $param_id=~/ACTIVE/) {
427	#The following pattern is used in several of these cases. It
428	#should probably be a subroutine...
429	#1. Get rid of double quotes in the data string
430	#2. Seperate the string into an array (each character is an entry)
431	#3. If it is empty put in a dash as the first character (not sure
432	# why I chose that...
433	#4. Fill the unused portion of the 50 character array with spaces
434	#5. Add the result to the master hash array for this param_id
435
436	#1
437	$value =~s/\"//g;
438	#2
439	@value = split("","$value");
440	@this_val = @value;
441	#3
442	unless ($#value+1) {
443	@value = ("-");
444	}
445	#4
446	for ($ii=$#value+1; $ii<50; $ii++) {
447	@this_val[$ii] = " ";
448	}
449	#5
450	push(@{$Vars{$param_id}},@this_val);
451
452	$expected_vars = $value;
453	$expected_vars=~s/\"//g;
454	@expected_vars = split('\s',$expected_vars);
455
456
457	my($t_type) = $param_id;
458	$t_type =~s/^ACTIVE\_//;
459
460	@{$expected_vars{$t_type}} = @expected_vars;
461	}
462	elsif ($type=~/CPU/) {
463	#see if ($type=~/SETTING/ && $param_id=~/ACTIVE/) {
464	$value =~s/\"//g;
465	@value = split("","$value");
466	@this_val = @value;
467	unless ($#value+1) {
468	@value = (" ");
469	}
470	for ($ii=$#value+1; $ii<50; $ii++) {
471	@this_val[$ii] = " ";
472	}
473	push(@{$Vars{$param_id}},@this_val);
474	}
475	elsif ($type=~/CLOCK/) {
476	#see if ($type=~/SETTING/ && $param_id=~/ACTIVE/) {
477	$value =~s/\"//g;
478	@value = split("","$value");
479	@this_val = @value;
480	unless ($#value+1) {
481	@value = (" ");
482	}
483	for ($ii=$#value+1; $ii<50; $ii++) {
484	@this_val[$ii] = " ";
485	}
486	push(@{$Vars{$param_id}},@this_val);
487	}
488
489	#The next little clause here is sort of weird. There are some
490	#computed or parsed variables that are not in the raw file as
491	#a param_id. For example PIRA PIRB and PIRC are all parameters
492	#but we really want to know what PIR is and we don't care about
493	#the raw data that goes into it... I'm just letting the script
494	#know to expect these variables to be included later.
495
496	elsif ($#{$expected_vars{$big_type}}) {
497	push(@{$loop_vars{$big_type}},$param_id);
498
499	if ($param_id=~/LINK$/) {
500	if ($first_one) {
501	@expected_vars = ('LINK','LINK_AMPS','LINK_VOLTS','LINK_AMP_HRS',
502	'LINK_TEMP');
503	push(@{$expected_vars{$big_type}},@expected_vars);
504
505	@expected_vars = ('SAMPLE_LOOP','BBAD_SAMPLE_LOOP');
506	push(@{$expected_vars{TIME}},@expected_vars);
507
508	@expected_vars = ('PIR');
509	push(@{$expected_vars{RLCAD}},@expected_vars);
510
511	$first_one = 0;
512	}
513	push(@{$loop_vars{$big_type}}, 'LINK','LINK_AMPS','LINK_VOLTS',
514	'LINK_AMP_HRS','LINK_TEMP');
515	}
516	}
517
518
519	#If we are at the start of a data record, use this time for the
520	#rest of the samples in this loop (may want to change this...)
521	if ($param_id=~/ATEMP/) {
522	$time_index++;
523	#Here I'm trying to work in some logic to deal with times when
524	#not every expected variable gets reported in the loop, for
525	#example when the BBAD has a time-out error. This leads to
526	#the system not reporting WSPDA or WDIRA, which is otherwise
527	#not noticed, and the empty values lead to a segmentation fault
528	#when the netcdf file is written.
529	foreach $expected_type (keys %expected_vars) {
530	$num_expected = $#{$expected_vars{$expected_type}};
531	$num_recieved = $#{$loop_vars{$expected_type}};
532
533	#Check to see if the number expected is greater than the number
534	#recieved (and also greater than 0, to deal with system info that
535	#doesn't describe variables in the loop!)
536	if ($num_expected>$num_recieved && $num_expected && $time_index>=0) {
537
538	#For every parameter that we are expecting, we should put in a
539	#fill value. Actually this all should probably go outside the
540	#loop as part of the initialization process.
541	foreach $expected_var (@{$expected_vars{$expected_type}}) {
542	#write fill value!
543
544	${$Vars{$expected_var}}[$time_index] = $fill_value;
545	} #End Foreach expected var
546
547	#Okay, this one is a bit wierd. Clock (sample_loop) is not
548	#yet included in RLCAD, BBAD, or SERIAL, so it isn't expected.
549	#Thus if it isn't reported, we need to insert a fill_value
550
551	#In the current version, this IS fixed. In fact CLOCK isn't used
552	#at all. It has been replaced by SAMPLE_LOOP.
553	unless(${$Vars{"CLOCK"}}[$time_index]) {
554	${$Vars{"CLOCK"}}[$time_index] = $fill_value;
555	} #end unless clock thing...
556	}#end if still expecting more...
557	}#end foreach expected type
558
559	%loop_vars = ();
560
561
562	#Now that is out of the way, time is a bit more normal, but still
563	#needs to be dealt with seperately as it occurs only once per
564	#sample loop and is not explicitly a loop variable.
565
566	$platform_time[$time_index] = $time;
567	$u_date = &UnixDate(&ParseDate("$time"),"%Y,%m,%d,%H,%M,%S");
568	($year,$month,$day,$hour,$minute,$second) = split(',',$u_date);
569
570	$u_date = &UnixDate(&ParseDate("$time UTC"),"%s");
571	$um1_date = &UnixDate(&ParseDate("$platform_time[$time_index -1] UTC"),
572	"%s");
573	$sample_interval = $u_date - $um1_date;
574	$hrs_per_sample = $sample_interval/3600; #Seconds/hr
575
576	unless ($time_index) {
577	$hrs_per_sample = 1/30;
578	}
579	#1970/01/01 00:00:00 as a datenum
580	$epoch_dn = 719529;
581
582	#Convert the seconds since epoch to a datenum
583	$datenum = $u_date/86400 + $epoch_dn;
584
585	${$Vars{'year'}}[$time_index] = $year;
586	${$Vars{'month'}}[$time_index] = $month;
587	${$Vars{'day'}}[$time_index] = $day;
588	${$Vars{'hour'}}[$time_index] = $hour;
589	${$Vars{'minute'}}[$time_index] = $minute;
590	${$Vars{'second'}}[$time_index] = $second;
591	${$Vars{'time'}}[$time_index] = $datenum;
592
593	}#end if ATEMP
594
595	#Apply information from the parse table to split up the output
596	&parse_vars;
597
598	if ($param_id=~/GPS1/) {
599	$param_id = 'GPSLON';
600	&convert_vars;
601	$param_id = 'GPSLAT';
602	&convert_vars;
603	$param_id = 'GPS1';
604	}
605
606
607	#Use the coefficients and formulas from the conversions table to
608	#apply any necessary conversions
609	&convert_vars;
610
611	#I'm not sure of a good way to get it to automatically convert
612	#and derive pir from the case temperature and the thermistor...
613	#This way is definitely a hack...
614
615	#If we have data for PIRA and PIRC, then we can compute PIR
616	if ($param_id =~/PIRC/ \| $param_id =~/PIRA/) {
617
618	if ($Vars{PIRA}[-1]!=~/$fill_value/ &&
619	$Vars{PIRC}[-1]!=~/$fill_value/){
620	$temporary_param = $param_id;
621	$param_id = 'PIR';
622
623	$pir_ct = $Vars{PIRA}[-1];
624	$pir_th = $Vars{PIRC}[-1];
625
626	&convert_vars;
627	$param_id = $temporary_param;
628	}
629	else {
630	$Vars{PIR} = $fill_value;
631	}
632	}
633	$value = undef;
634
635	}
636
637	#Otherwise the data line has no commas, so it must be an adcp line
638	else {
639	#I haven't put any logic in here yet
640	#&parse_adcp;
641	}
642	}
643
644	######################################################################
645	### Parse Vars #######################################################
646	sub parse_vars {
647
648	#Seperate out the output from matching conditions
649
650	#The hash array parse_output comes from the meta_data perl module that
651	#looks into the data tables. It explains the string format for each
652	#variable and how to assign the variables given that format.
653	$output = $parse_output{$param_id};
654	@output = split(',',$output);
655
656
657	if ($units{$param_id}=~/string/) {
658	#It was complaining about the gps checksum until I added ''
659	eval("$output[0] = '$value';");
660	}
661	else {
662	#Apply matching conditions from the parse table
663	$old_v = $value;
664
665	#Okay, this is hard wired... I know it's better to avoid, but I wanted
666	#a quick fix... This handles wind, which is an array of values separated
667	#by spaces. An average value should be presented for later conversion.
668	if ($parse_info{$param_id}=~/\@avg/) {
669	#Spit up space delimited values and ignore fill values
670	$old_v =~s/\ \-9999//g;
671	@value = split(' ',$old_v);
672	eval("$output[0]".' = @value;');
673	}
674	elsif ($parse_info{$param_id}=~/\@gva/) {
675	#This is for the compass pitch and roll which is weird...
676	@comp = @pitch = @roll = ();
677
678	@value = split('\$',$old_v);
679	foreach $vals (@value) {
680	$vals =~m/$C([0-9\.\-])P([0-9\.\-])R([0-9\.\-])T./;
681	push(@comp,$1);
682	push(@pitch,$2);
683	push(@roll,$3);
684	#print "so $1, $2, $3\n";
685	}
686	#eval("$output[0]".' = @value;');
687	}
688
689	#If the output is one parameter (applies to almost all sensors)
690	elsif ($#output==0) {
691	#Assign the output value to the parse_info matching condition
692	#Here $1 signifies the part that matches what is inside the
693	#first set of parenthesis in a matching condition.
694	#Example: if $value is "PS=+14.693"
695	# then $value =~m/"PS\=\+(.*)"/;
696	# sets $1 to be 14.693
697	$value =~m/$parse_info{$param_id}/;
698	@value = $1;
699
700	#As long as there is a defined value, set the output to that
701	if (defined $value[0]) {
702	eval("$output[0]".' = $value[0];');
703	}
704	else {
705	eval("$output[0] = $fill_value;");
706	}
707
708	}
709
710	#Applies to the ctd and the link, in particular
711	elsif ($#output>0) {
712	if ($param_id =~m/CTD/) {
713	#As long as the string begins with '"#0' proceed.
714	if ($value=~m/^\"\#0/) {
715	$pi = $parse_info{$param_id};
716	}
717	else {
718	$pi = 'NOBAD!!';
719	}
720	}
721	else {
722
723	$pi = $parse_info{$param_id};
724	}
725	if ($param_id=~/LINK/) {
726	#print "$value\n$pi\n";
727	}
728	elsif ($param_id=~/GPS1/) {
729	#print "$value\n";
730	@v = split('',$value);
731	$mygps_lat = join('',@v[11..17]);
732	$mygps_lon = join('',@v[19..25]);
733
734	}
735	$value =~m/$pi/;
736
737	#Loop through each output parameter and assign the value
738	#For example suppose there were three variables in the matching
739	#condition '$value =~m/$pi/;', then they will be assigned to
740	#$1,$2, and $3, thus you want to assign them as
741	#eval("$output[0] = $1;");
742	#so if @output = ($temp,$humidity,$pressure)
743	#then you have $temp set to the contents of $1. Clear???
744	for ($out_num=1; $out_num<=$#output+1; $out_num++) {
745	eval("$output[$out_num-1] =\$$out_num;");
746	}
747	if ($param_id=~/LINK/) {
748	#print ": $link_time, $link_volts, \n";
749	}
750	elsif ($param_id=~/GPS1/) {
751	if ($mygps_lat+0) {$gps_lat = $mygps_lat;}
752	else {$gps_lat = $fill_value;}
753	if ($mygps_lon+0) {$gps_lon = $mygps_lon;}
754	else {$gps_lon = $fill_value;}
755	}
756	}
757	}
758	}
759
760	######################################################################
761	### Convert Vars #####################################################
762	sub convert_vars {
763	our($fill_value);
764
765	#Again, these coefficients and formulas are coming from the
766	#meta_data perl module, and the corresponding data tables.
767
768	#Since conversion coefficients are different for the different
769	#platforms it looks like this $R4_conversion_coeff{PSP}
770
771	#Get conversion coefficients and formula from hash array for
772	#each tower id and each parameter id.
773	eval ('$c_coeff = $'.$tower_id.'_conversion_coeff{$param_id};');
774	eval ('$c_formula = $'.$tower_id.'_conversion_formula{$param_id};');
775
776	#Split up conversion coefficients from the converstions table
777	@v = split(',',$c_coeff);
778
779	#If there is any set of conversion coefficients (and hence a
780	#conversion to be applied), apply them!
781	if ($#v+1) {
782	$temp_var = '';
783	if ($c_formula=~/^\&/) {
784	#If the conversion formula is a subroutine call, evaluate!
785	eval("$c_formula");
786	}
787	else {
788	#As long as the output wasn't the fill value, convert it.
789	my ($out) = eval("$output[0]");
790	unless($out=~/$fill_value/) {
791	$temp_var = eval("$c_formula");
792	}
793
794	#Otherwise, we don't want to convert a fill value to
795	#something else!
796	else {
797	$temp_var = $fill_value;
798	}
799
800	#Set the array index based upon the record dimension
801	&set_index;
802
803	#This should go somewhere else...
804	if ($param_id=~/RAIN/) {
805	#print "Hrs per : $hrs_per_sample\n";
806	$temp_var = $temp_var/$hrs_per_sample;
807	}
808
809	#Finally, set this value of the Vars super duper array
810	#for this time to the value that we've computed
811	unless ($index +1) {
812	print "Warning. $param_id has index $index. Adjusting...\n";
813	${$Vars{$param_id}}[$index+1] = $temp_var;
814	}
815	else {
816	${$Vars{$param_id}}[$index] = $temp_var;
817	}
818	}
819	}
820
821	else {
822	#print "No conversion info available for $param_id at $tower_id\n";
823	}
824	}
825
826
827	######################################################################
828	### Handle Input #####################################################
829	sub handle_input {
830
831	#By putting default values in a file, the script should easily switch
832	#back and forth between the afs version and the local nemo version!
833	open(DEFAULTS,"${bin_dir}tower_parse.defaults");
834	while (<DEFAULTS>) {
835	#Get rid of the newline
836	chop;
837	#Ignore comments starting with #
838	$_=~s/#.*//;
839	#If the line is non-empty
840	if ($_) {
841	my($var_name, $default_value) = split(',');
842	eval("$var_name = $default_value;");
843	}
844	}
845
846	($date_string,%varargin) = @ARGV;
847
848	#Associate each input argument pair with the appropriate variable
849	%varargout = (
850	'-out' => '$out_file_name',
851	'-d' => '$debug',
852	'-t' => '$tower_id',
853	'-p' => '$package',
854	);
855
856	#Set each variable to the value of the associated input argument
857	#Example 1: '-d 1' sets the debug flag to true
858	#Example 2: '-out 2004_06.nc' sets the output file name to 2004_06.nc
859	foreach $argin (keys %varargin) {
860	eval("$varargout{$argin} = '$varargin{$argin}';");
861	}
862
863	$d_data_dir = "$d_data_dir$tower_id/";
864	$out_data_dir = "$out_data_dir$tower_id/";
865
866	#Since there are a few different types of date string inputs, we need
867	#to figure out how to deal with this one!
868	if ($date_string) {
869	&handle_date;
870	}
871
872	#If there's no date string, then we don't know what to look for!
873	else {
874	die "Must at least input a date string!\n".
875	"./tower_parse.perl '200406' -out 2004_06.nc -d 1\n";
876	}
877
878	#Helps the user chooose an appropriate filename
879	unless ($out_file_name=~/.nc$/) {
880	if ($out_file_name) {
881	$out_file_sname = "${out_file_name}_sys.nc";
882	$out_file_name = "${out_file_name}.nc";
883	}
884	else {
885	$out_file_name = $d_out_file_name;
886	@do = split('',$d_out_file_name);
887	$out_file_sname = join('',@do[0..$#do-3],'_sys.nc');
888	}
889	}
890	else {
891	unless ($out_file_sname) {
892	@do = split('',$out_file_name);
893	$out_file_sname = join('',@do[0..$#do-3],'_sys.nc');
894	}
895	}
896
897	#This will come up many times throughout the script. Another way of
898	#looking at the debug flag is that it turns on comments (verbose).
899	if ($debug) {
900	print "\n\nData dir : $d_data_dir\nOutput dir : $out_data_dir\n".
901	"Outfile name : $out_file_name and $out_file_sname\n";
902	if ($#dates) {
903	print "\nLooking for data between $dates[0] and $dates[-1]\n";
904	}
905	}
906
907
908	}
909
910
911	######################################################################
912	### Handle Date ######################################################
913	sub handle_date {
914	#Deal with date input.
915	if ($date_string=~/\ /) {
916	($st_date,$en_date) = split(' ',$date_string);
917	@st_date = split('',$st_date);
918	@en_date = split('',$en_date);
919
920	#print "So $st_date, $en_date is $#st_date or @st_date\n";
921	$st_date = join('',(@st_date[0..3],' ',@st_date[4..5],' ',@st_date[6..7]));
922	$en_date = join('',(@en_date[0..3],' ',@en_date[4..5],' ',@en_date[6..7]));
923
924	#print "So $st_date, $en_date is @st_date\n";
925	$this_date = $st_date;
926	$this_parsed = &ParseDate($st_date);
927	$en_parsed = &ParseDate($en_date);
928
929	@dates = ();
930	%month_dirs = ();
931
932	while ($this_parsed <= $en_parsed) {
933	$ymd = &UnixDate($this_parsed,"%Y%m%d");
934	push(@dates,$ymd);
935	$month_dirs{$ymd} = &UnixDate($this_parsed,"%Y_%m");
936	$this_date = &UnixDate(&DateCalc($this_date,"+ 1 day"), "%Y %m %d");
937	$this_parsed = &ParseDate($this_date);
938	}
939	#print "$#dates Dates! @dates\n";
940	}
941	elsif ($date_string=~/latest/) {
942	$today = &ParseDate("today");
943	$date1 = &UnixDate(&DateCalc("$today","- 2 day"),"%Y%m%d");
944	$date2 = &UnixDate(&DateCalc("$today","- 1 day"),"%Y%m%d");
945	$date3 = &UnixDate("$today","%Y%m%d");
946	$mond1 = &UnixDate(&DateCalc("$today","- 2 day"),"%Y_%m");
947	$mond2 = &UnixDate(&DateCalc("$today","- 1 day"),"%Y_%m");
948	$mond3 = &UnixDate("$today","%Y_%m");
949	@dates = ($date1,$date2,$date3);
950	%month_dirs = ($date1,$mond1,$date2,$mond2,$date3,$mond3);
951	#print "@dates\n";
952	}
953	elsif ($date_string =~/help/) {
954	#oops this isn't a date string at all!
955	# I should probably put this somewhere else...
956	print "\nUsage: ./tower_parse.perl date_string\n\n".
957	"Optional argument pairs:\n".
958	" -out output file name\n".
959	" -d debug\n".
960	# " -dir directory path to data\n".
961	" -t tower id\n".
962	" -nemo scp data file to nemo\n\n";
963	print "Date string can have any of the following formats:\n".
964	" yyyymm Finds all data for the month\n".
965	" yyyymmdd Finds data just for one day\n".
966	" 'yyyymmdd yyyymmdd' Finds data within the date range\n".
967	" latest Finds today, yesterday and the previous day\n";
968	die "\n";
969	}
970	else {
971	@dates = ($date_string);
972	@date_str = split('',$date_string);
973	$fn = join('',(@date_str[0..3],'_',@date_str[4..5]));
974
975	unless($out_file_name) {
976	#Unless it has been input by the user!
977	$out_file_name = "${tower_id}_${fn}.nc";
978	$out_file_sname = "${tower_id}_${fn}_sys.nc";
979	}
980
981	#Also change default data directory!
982	$month_dirs{$date_string} = $fn;
983	}
984	@month_dirs = values (%month_dirs);
985	}
986
987
988	######################################################################
989	### Parse ADCP #######################################################
990	sub parse_adcp {
991
992	#print "We have ADCP data!\n";
993	#$data = join('',@data);
994
995	}
996
997
998	######################################################################
999	### Vardef Attput ####################################################
1000	sub vardef_attput {
1001
1002	#This is the generic variable definition and attribute put used by
1003	#write_netcdf.
1004	$format = 'NetCDF::'.$format{$param};
1005	if ($debug>2) {
1006	print "$param is of format $format\n";
1007	}
1008
1009	@dimensions = split('\s',$dimensions{$param});
1010	@dim_names = @dim_names{@dimensions};
1011	@dim_ids = @dimid{@dim_names};
1012
1013	$dimensions = join('',@dimensions);
1014
1015	#Set the record dimension to be time or file number based on
1016	#contents of variables.table
1017	if ($dimensions=~/f/) {
1018	$ncid{$param} = $ncid_system;
1019	}
1020	elsif ($dimensions=~/t/) {
1021	$ncid{$param} = $ncid;
1022	}
1023
1024
1025	$varid{$param} = NetCDF::vardef($ncid{$param},$param,
1026	eval("$format"),
1027	[@dim_ids]);
1028
1029	die "Couldn't define variable $param\n" if $varid{$param}<0;
1030
1031
1032	#print "$param should have $#all_attributes\n";
1033	foreach $attribute (@all_attributes) {
1034	$att_val = eval('\$'.$attribute.'{$param}');
1035	if (defined $att_val && $att_val!=~/^N\/A/) {
1036	$attid = NetCDF::attput($ncid{$param},$varid{$param},"$attribute",
1037	NetCDF::CHAR,
1038	eval('\$'.$attribute.'{$param}'));
1039	}#end if
1040	}#end foreach $attribute
1041	}
1042
1043	######################################################################
1044	### Varput ###########################################################
1045	sub varput {
1046	$format = 'NetCDF::'.$format{$param};
1047
1048	@dimensions = split('\s',$dimensions{$param});
1049	@dim_names = @dim_names{@dimensions};
1050	#@dim_ids = @dimid{@dim_names};
1051	@dim_counts = @dims{@dim_names};
1052	@dim_starts = @zeros[0..$#dim_names];
1053
1054	unless ($param =~/CLOCK/ \|\| $param =~/BBAD/ \|\|
1055	$param =~/GPS$/ \|\| $param =~/LINK$/) {
1056
1057	#No matter how good it looks, test the arrray for undefined values!
1058	#(This is painfully slow and really should be replaced by initializing
1059	#the arrays at the start...)
1060	@temp_array = @{$Vars{$param}};
1061	for ($ii=0; $ii<=$#temp_array; $ii++) {
1062	unless (defined($temp_array[$ii])) {
1063	${$Vars{$param}}[$ii] = $fill_value;
1064	}
1065	}
1066
1067	if ($debug>1) {
1068	print "$param length ".eval($#{$Vars{$param}}+1).
1069	", id $varid{$param}\n";
1070	}
1071
1072	unless ($format=~/CHAR/) {
1073	#Put the variable into the netcdf file!
1074	$status = NetCDF::varput($ncid{$param},$varid{$param},
1075	\@dim_starts,\@dim_counts,
1076	\@{$Vars{$param}});
1077	}
1078	else {
1079	#Test to make sure it is as long as we expected...
1080	if ($#{$Vars{$param}}+1==($#tower_files+1)*50) {
1081	$status = NetCDF::varput($ncid{$param},$varid{$param},
1082	\@dim_starts,\@dim_counts,
1083	\@{$Vars{$param}});
1084	}
1085	elsif ($debug>1) {
1086	print "Didn't have enough values for $param\n";
1087	}
1088	}
1089	}#End unless that ignores clock,bbad,gps and link strings
1090	else {
1091	print "$param being ignored\n";
1092	}
1093	}
1094
1095	######################################################################
1096	### Compute Sub Vars #################################################
1097	sub compute_sub_vars {
1098
1099	if ($#sub_vars +1) {
1100	$orig_param = $param_id;
1101
1102	foreach $sub_param (@sub_vars) {
1103	$param_id = $sub_param;
1104	&convert_vars;
1105	}
1106	$param_id = $orig_param;
1107	}
1108	}
1109
1110	######################################################################
1111	### Set Index ########################################################
1112	sub set_index {
1113
1114	#If the record dimension is time, use the time index
1115	if ($dimensions{$param_id}=~/t/) {
1116	$index = $time_index;
1117	}
1118	#If the record dimension is file number use that index
1119	elsif ($dimensions{$param_id}=~/f/) {
1120	$index = $file_index;
1121	}
1122	}
1123
1124
1125	######################################################################
1126	### Avg Wind #########################################################
1127	sub avg_wind {
1128	#1. Check for all necessary components:
1129	# a. 'spda' -> WSPDA
1130	# b. 'dira' -> WDIRA
1131	# c. 'spdb' -> WSPDB (Optional) [alternative is 'spdn' if empty]
1132	# d. 'dirb' -> WDIRB (Optional) [alternative is 'dirn' if empty]
1133	# e. 'cmp1' -> COMP1 (Optional) [alternative is an offset, which
1134	# defaults to 0 unless otherwise set... Does not requre pkgb]
1135	# f. 'cmp2' -> COMP2 [Not Used!]
1136	#
1137	#2. Loop through times:
1138	# a. Replace out of range values or empty package with fill_value
1139	# b. Apply direction offset
1140	# c. Compute vector components (u,v)
1141	#
1142	#3. Compute statistics
1143	#
1144	#4. Store data values
1145
1146	my(@junk) = ();
1147	my(@input) = (@_);
1148	%names = ('spda'=>'wi_spd_a',
1149	'dira'=>'wi_dir_a',
1150	'spdb'=>'wi_spd_b',
1151	'dirb'=>'wi_dir_b',
1152	'spdn'=>'junk',
1153	'dirn'=>'junk',
1154	'cmp1'=>'comp_1',
1155	'cmp2'=>'comp_2');
1156
1157
1158	#1. Check for all necessary components:
1159	if ($input[0]=~/spda/) {
1160	#First variable on the list, initialize temporary data arrays
1161
1162	#If we've already been through the loop once...
1163	if ($dira[0]) {
1164	$real_last_var = $last_var;
1165
1166	if ($#cmp1<0) {
1167	#No compass data. Nor is there an available offset.
1168	@cmp1 = (0);
1169	}
1170	}
1171	else {
1172	#for the first one, so it doesn't confuse the real_last_var
1173	$real_last_var = 'NOPE';
1174	}
1175
1176	#Now initialize the arrays
1177	@spda = @dira = @spdb = @dirb = @cmp1 = @cmp2 = ();
1178	}
1179	elsif ($input[0]=~/dir/) {
1180	if ($#input==1) {
1181	#Then we have been given a directional offset and should
1182	#apply it!
1183	$offset = $input[1];
1184	}
1185	else {
1186	$offset = 0;
1187	}
1188	}
1189	#print("\@$input[0] = \@$names{$input[0]};\n");
1190	eval("\@$input[0] = \@$names{$input[0]};");
1191
1192	if ($input[0]=~/$real_last_var/) {
1193	#Write the data from this sample.
1194	&handle_wind;
1195	&write_wind;
1196	}
1197
1198	$last_var = $input[0];
1199	}
1200
1201
1202	######################################################################
1203	### Handle Wind ######################################################
1204	sub handle_wind {
1205	#This subroutine works for average_wind: once data arrays are
1206	#collected it preforms all the necessary operations on them.
1207
1208	my($time);
1209	our($u_mean_a,$u_max_a,$u_std_a,$u_mean_b,$u_max_b,$u_std_b,
1210	$v_mean_a,$v_max_a,$v_std_a,$v_mean_b,$v_max_b,$v_std_b,
1211	$mean_spd_a,$mean_dir_a,$mean_spd_b,$mean_dir_b,$spd_std_a,
1212	$dir_std_a,$spd_std_b,$dir_std_b,$gust_a,$gust_b) = $fill_value;
1213
1214	my(@u_a,@v_a,@u_b,@v_b,@temp_cmp,@comp_off) = ();
1215
1216
1217	#Loop through samples...
1218	for ($time=0; $time<=$#spda; $time++) {
1219
1220	#Make sure there is some directional offset to apply!
1221	if ($#cmp1<1) {
1222	#If compass data wasn't parsed correctly, fix it!
1223	if ($cmp1[0]=~/\ /) {
1224	@temp_cmp = split(' ',$cmp1[0]);
1225	@cmp1 = @temp_cmp;
1226	$comp_off[$time] = &volts2comp($cmp1[$time]) + $offset;
1227	}
1228	#No compass, using fixed offset
1229	else {
1230	$comp_off[$time] = &volts2comp($cmp1[0]) + $offset;
1231	}
1232	}
1233	else{
1234	#There is a compass reading so use that (buoy)
1235	$comp_off[$time] = &volts2comp($cmp1[$time]) + $offset;
1236	}
1237
1238
1239	#Apply simple, static conversion from voltage to scientific values
1240	($spda[$time],$dira[$time]) = &volts2spddir($spda[$time],$dira[$time]);
1241	($spdb[$time],$dirb[$time]) = &volts2spddir($spdb[$time],$dirb[$time]);
1242
1243
1244	#Make sure speeds are reasonable and if speed or direction are
1245	#bad, both are badflagged
1246	$bada = ($spda[$time]>100 \| $spda[$time]==$fill_value \|
1247	$dira[$time]==$fill_value);
1248	if ($bada) {
1249	$spda[$time] = $dira[$time] = $u_a[$time] = $v_a[$time] = $fill_value;
1250	}
1251	else {
1252	#Take the direction mod 360
1253	$dira[$time] = &modulo($dira[$time] + $comp_off[$time], 360);
1254	#Convert Spd and Dir to U and V components
1255	($u_a[$time],$v_a[$time]) = &spddir2uv($spda[$time],$dira[$time]);
1256	}
1257
1258	#Ditto for the b package! (as long as there is a b package...)
1259	unless ($#spdb<0) {
1260	$badb = ($spdb[$time]>100 \| $spdb[$time]==$fill_value \|
1261	$dirb[$time]==$fill_value);
1262	if ($badb) {
1263	$spdb[$time] = $dirb[$time] = $u_b[$time] =
1264	$v_b[$time] = $fill_value;
1265	}
1266	else {
1267	$dirb[$time] = &modulo($dirb[$time] + $comp_off[$time], 360);
1268	($u_b[$time],$v_b[$time]) = &spddir2uv($spdb[$time],$dirb[$time]);
1269	}
1270	}
1271	}
1272
1273	#a package stats
1274	#($mean,$std,$max) = &fill_stats(@array);
1275	($u_mean_a,$u_std_a,$u_max_a) = &fill_stats(@u_a);
1276	($v_mean_a,$v_std_a,$v_max_a) = &fill_stats(@v_a);
1277	#In terms of speed, mean speed is wrong
1278	(undef,$spd_std_a,$gust_a) = &fill_stats(@spda);
1279	#For direction, only std dev is right...
1280	(undef,$dir_std_a,undef) = &fill_stats(@dira);
1281
1282	($mean_spd_a,$mean_dir_a) = &uv2spddir($u_mean_a,$v_mean_a);
1283
1284	($comp_mean,$comp_std,undef) = &fill_stats(@comp_off);
1285
1286	#b package stats
1287	unless ($#u_b<0) {
1288	($u_mean_b,$u_std_b,$u_max_b) = &fill_stats(@u_b);
1289	($v_mean_b,$v_std_b,$v_max_b) = &fill_stats(@v_b);
1290	#In terms of speed, mean speed is wrong
1291	(undef,$spd_std_b,$gust_b) = &fill_stats(@spdb);
1292	#For direction, only std dev is right...
1293	(undef,$dir_std_b,undef) = &fill_stats(@dirb);
1294
1295	($mean_spd_b,$mean_dir_b) = &uv2spddir($u_mean_b,$v_mean_b);
1296	}
1297	}
1298
1299	######################################################################
1300	### Write_Wind ########################################################
1301	sub write_wind {
1302
1303	unless ($mean_dir_a == $fill_value \|\|
1304	$mean_spd_a == $fill_value) {
1305	#put variables into array!!!
1306	${$Vars{WSPDA}}[$index] = $mean_spd_a;
1307	${$Vars{WDIRA}}[$index] = $mean_dir_a;
1308	${$Vars{WEASTA}}[$index] = $u_mean_a;
1309	${$Vars{WNORTHA}}[$index] = $v_mean_a;
1310	${$Vars{WSPDSTDA}}[$index] = $spd_std_a;
1311	${$Vars{WDIRSTDA}}[$index] = $dir_std_a;
1312	${$Vars{WGSTA}}[$index] = $gust_a;
1313	}
1314	else {
1315	${$Vars{WSPDA}}[$index] = $fill_value;
1316	${$Vars{WDIRA}}[$index] = $fill_value;
1317	${$Vars{WEASTA}}[$index] = $fill_value;
1318	${$Vars{WNORTHA}}[$index] = $fill_value;
1319	${$Vars{WSPDSTDA}}[$index] = $fill_value;
1320	${$Vars{WDIRSTDA}}[$index] = $fill_value;
1321	${$Vars{WGSTA}}[$index] = $fill_value;
1322	}
1323	unless ($mean_dir_b == $fill_value \|\|
1324	$mean_spd_b == $fill_value) {
1325	#put variables into array!!!
1326	${$Vars{WSPDB}}[$index] = $mean_spd_b;
1327	${$Vars{WDIRB}}[$index] = $mean_dir_b;
1328	${$Vars{WEASTB}}[$index] = $u_mean_b;
1329	${$Vars{WNORTHB}}[$index] = $v_mean_b;
1330	${$Vars{WSPDSTDB}}[$index] = $spd_std_b;
1331	${$Vars{WDIRSTDB}}[$index] = $dir_std_b;
1332	${$Vars{WGSTB}}[$index] = $gust_b;
1333	}
1334	else {
1335	${$Vars{WSPDB}}[$index] = $fill_value;
1336	${$Vars{WDIRB}}[$index] = $fill_value;
1337	${$Vars{WEASTB}}[$index] = $fill_value;
1338	${$Vars{WNORTHB}}[$index] = $fill_value;
1339	${$Vars{WSPDSTDB}}[$index] = $fill_value;
1340	${$Vars{WDIRSTDB}}[$index] = $fill_value;
1341	${$Vars{WGSTB}}[$index] = $fill_value;
1342	}
1343
1344	unless ($comp_std == 0 \|\|
1345	$comp_mean == -9999) {
1346	${$Vars{COMP1}}[$index] = $comp_mean;
1347	}
1348	else {
1349	${$Vars{COMP1}}[$index] = $fill_value;
1350	}
1351	}
1352
1353	######################################################################
1354	### Handle Comp ######################################################
1355	sub handle_comp {
1356	#This subroutine works for average_wind: once data arrays are
1357	#collected it preforms all the necessary operations on them.
1358
1359	#Trying to figure out what to do with hourly compass measurements
1360	#including pitch and roll. Need mean and std
1361
1362	#print "HANDLE COMP!$#comp $#pitch $#roll\n";
1363
1364	our($compd_mean,$compd_std,$compd_max,$pitch_mean,$pitch_std,
1365	$pitch_max,$roll_mean,$roll_std,$roll_max) = $fill_value;
1366
1367	my(@u_c,@v_c,@u_p,@v_p,@u_r,@v_r) = ();
1368
1369	#print "Okay? $comp[0] $pitch[0] $roll[0]?\n";
1370	if ($#comp>0) {
1371	#print "ARRAY: @comp\n";
1372	#print "Looping compass...\n";
1373	#Loop through samples...
1374	for ($time=0; $time<=$#comp; $time++) {
1375	#Since we don't have a speed, we will use the unit vector for
1376	#the averaging...
1377
1378	#print "C: $comp[$time] P: $pitch[$time] R: $roll[$time] at 1\n";
1379
1380	#Take the direction mod 360
1381	$comp[$time] = &modulo($comp[$time], 360);
1382	$pitch[$time] = &modulo($pitch[$time], 360);
1383	$roll[$time] = &modulo($roll[$time], 360);
1384
1385	#print "C: $comp[$time] P: $pitch[$time] R: $roll[$time] at 2\n";
1386
1387	#Convert Spd and Dir to U and V components
1388	($u_c[$time],$v_c[$time]) = &spddir2uv(1,$comp[$time]);
1389	($u_p[$time],$v_p[$time]) = &spddir2uv(1,$pitch[$time]);
1390	($u_r[$time],$v_r[$time]) = &spddir2uv(1,$roll[$time]);
1391
1392	#print "C: $u_c[$time] $v_c[$time] P: $u_p[$time] $v_p[$time]".
1393	#" R: $u_r[$time] $v_r[$time] at 3\n";
1394	}
1395
1396	#($mean,$std,$max) = &fill_stats(@array);
1397	($u_c_mean,$u_c_std,$u_c_max) = &fill_stats(@u_c);
1398	($u_p_mean,$u_p_std,$u_p_max) = &fill_stats(@u_p);
1399	($u_r_mean,$u_r_std,$u_r_max) = &fill_stats(@u_r);
1400
1401	($v_c_mean,$v_c_std,$v_c_max) = &fill_stats(@v_c);
1402	($v_p_mean,$v_p_std,$v_p_max) = &fill_stats(@v_p);
1403	($v_r_mean,$v_r_std,$v_r_max) = &fill_stats(@v_r);
1404
1405	#Now get back to the direction!
1406	(undef,$compd_mean) = &uv2spddir($u_c_mean,$v_c_mean);
1407	(undef,$compd_std) = ($u_c_std2+$v_c_std2)**(1/2);
1408	(undef,$compd_max) = &uv2spddir($u_c_max,$v_c_max);
1409
1410	(undef,$pitch_mean) = &uv2spddir($u_p_mean,$v_p_mean);
1411	(undef,$pitch_std) = ($u_p_std2+$v_p_std2)**(1/2);
1412	(undef,$pitch_max) = &uv2spddir($u_p_max,$v_p_max);
1413
1414	(undef,$roll_mean) = &uv2spddir($u_r_mean,$v_r_mean);
1415	(undef,$roll_std) = ($u_r_std2+$v_r_std2)**(1/2);
1416	(undef,$roll_max) = &uv2spddir($u_r_max,$v_r_max);
1417
1418
1419	#print "Compass: $comp_mean, $comp_std, $comp_max\n@comp\n";
1420	#print "Pitch: $pitch_mean, $pitch_std, $pitch_max\n@pitch\n";
1421	#print "Roll: $roll_mean, $roll_std, $roll_max\n@roll\n";
1422
1423	${$Vars{COMPMEAN}}[$file_index] = $compd_mean;
1424	${$Vars{COMPSTD}}[$file_index] = $compd_std;
1425	${$Vars{COMPMAX}}[$file_index] = $compd_max;
1426	${$Vars{PITCHMEAN}}[$file_index] = $pitch_mean;
1427	${$Vars{PITCHSTD}}[$file_index] = $pitch_std;
1428	${$Vars{PITCHMAX}}[$file_index] = $pitch_max;
1429	${$Vars{ROLLMEAN}}[$file_index] = $roll_mean;
1430	${$Vars{ROLLSTD}}[$file_index] = $roll_std;
1431	${$Vars{ROLLMAX}}[$file_index] = $roll_max;
1432	}
1433	else {
1434	#print "Nope...\n";
1435	${$Vars{COMPMEAN}}[$file_index] = $fill_value;
1436	${$Vars{COMPSTD}}[$file_index] = $fill_value;
1437	${$Vars{COMPMAX}}[$file_index] = $fill_value;
1438	${$Vars{PITCHMEAN}}[$file_index] = $fill_value;
1439	${$Vars{PITCHSTD}}[$file_index] = $fill_value;
1440	${$Vars{PITCHMAX}}[$file_index] = $fill_value;
1441	${$Vars{ROLLMEAN}}[$file_index] = $fill_value;
1442	${$Vars{ROLLSTD}}[$file_index] = $fill_value;
1443	${$Vars{ROLLMAX}}[$file_index] = $fill_value;
1444	}
1445	}
1446
1447
1448
1449	######################################################################
1450	### Volts2SpdDir #####################################################
1451	sub volts2spddir {
1452
1453	my($volts_spd,$volts_dir) = @_;
1454	my($wind_spd,$wind_dir) = undef;
1455
1456	unless($volts_spd=~/$fill_value/) {
1457	#0->1 V : 0 60 m/s
1458	$wind_spd = $volts_spd*60;
1459	}
1460	else { $wind_spd = $fill_value; }
1461
1462	unless($volts_dir=~/$fill_value/) {
1463	#0->1 V : 0 360 deg
1464	$wind_dir = $volts_dir*360;
1465	}
1466	else { $wind_dir = $fill_value; }
1467
1468	return ($wind_spd,$wind_dir);
1469	}
1470
1471	######################################################################
1472	### Volts2Comp #######################################################
1473	sub volts2comp {
1474
1475	my($volts_comp) = @_;
1476
1477	unless($volts_comp=~/$fill_value/) {
1478	#0->1 V : 0 360 deg
1479	return $volts_comp*360;
1480	}
1481	else { return $fill_value; }
1482	}
1483
1484
1485
1486	######################################################################
1487	### SpdDir2UV ########################################################
1488	sub spddir2uv {
1489	#convert speed and direction to u and v
1490	my($wind_spd,$wind_dir) = @_;
1491	my($deg2rad) = .0174533;
1492	my($wind_u,$wind_v) = undef;
1493
1494
1495	unless($wind_spd=~/$fill_value/ \|\| $wind_dir=~/$fill_value/) {
1496	#okay this is terrible but I'm at a loss what else to do.
1497	$wind_v = ($wind_spd)cos($wind_dir$deg2rad);
1498	$wind_u = ($wind_spd)sin($wind_dir$deg2rad);
1499	return ($wind_u,$wind_v);
1500	}
1501	else {
1502	#print "no!?";
1503	return ($fill_value,$fill_value);
1504	}
1505	}
1506
1507	######################################################################
1508	### UV2SpdDir ########################################################
1509	sub uv2spddir {
1510	#convert speed and direction to u and v
1511	my($wind_u,$wind_v) = @_;
1512	my($deg2rad) = .0174533;
1513	my($wind_spd,$wind_dir) = undef;
1514
1515	unless($wind_u=~/$fill_value/ \|\| $wind_v=~/$fill_value/) {
1516	$wind_spd = ($wind_u2 + $wind_v2)**(1/2);
1517	$wind_dir = &modulo(atan2($wind_u, $wind_v)/$deg2rad,360);
1518	return ($wind_spd,$wind_dir);
1519	}
1520	else {
1521	return ($fill_value,$fill_value);
1522	}
1523	}
1524
1525	######################################################################
1526	### Avg Rain #########################################################
1527	sub avg_rain {
1528	#Convert cumulative rainfall readings for self siphoning rain
1529	#guage into a rainfall amount over the sampling period.
1530
1531	#I haven't figured this out yet...
1532	#print "RAIN???\n";
1533
1534	#First I need to check the voltage -> volume.
1535	#Then, make sure it isn't in the process of dumping (look at diff,
1536	#confirm that it isn't largely negative)
1537	#Finally look at the change since the previous sample loop to get
1538	#the cumulative rain (and if it is negative, add the dumping
1539	#volume)
1540
1541	#print "Testing rain:\n";
1542
1543	#my ($ind,$raindiff);
1544	#for($ind=1;$ind<$#cumrain;$ind++) {
1545	# $raindiff = $cumrain[$ind]-$cumrain[$ind-1];
1546	# if ($raindiff<-.05) {
1547	# print "Dumping???\n";
1548	# }
1549	# elsif ($raindiff>.05) {
1550	# print "Raining???\n";
1551	# }
1552	#}
1553	my($maxrain);
1554	$maxrain = &max(@cumrain);
1555
1556
1557	unless ($maxrain=~/$fill_value/) {
1558	$maxrain = $maxrain*10;
1559	}
1560	else {
1561	$maxrain = $fill_value;
1562	}
1563	${$Vars{RAIN}}[$index] = $maxrain;
1564	}
1565
1566
1567
1568
1569	######################################################################
1570	### Fill Stats #######################################################
1571	sub fill_stats(@array) {
1572
1573	my(@array,$mean_value);
1574	@array = @_;
1575
1576	$array = join(',',@array);
1577	$array =~s/$fill_value//g;
1578	$array =~s/\,\,/\,/g;
1579	@array = split(',',$array);
1580
1581	unless ($#array<0) {
1582	($mean) = &mean(@array);
1583	($std) = &std_dev(@array);
1584	($max) = &max(@array);
1585	return ($mean,$std,$max);
1586	}
1587	else {
1588	return ($fill_value,$fill_value,$fill_value);
1589	}
1590	}
1591
1592	######################################################################
1593	### Mean #############################################################
1594	sub mean(@array) {
1595	my(@array,$mean_value);
1596	@array = @_;
1597
1598	$sum_array = join('+',@array);
1599
1600	if ($sum_array =~m/\+\+/) {
1601	return $fill_value;
1602	}
1603	#print "Averaging!\n$sum_array\n";
1604
1605	eval('$mean_value = ('."$sum_array".')/($#array+1);');
1606	return $mean_value;
1607	}
1608
1609	######################################################################
1610	### Std Dev ##########################################################
1611	sub std_dev(@array) {
1612	my(@array,$std_dev);
1613	@array = @_;
1614
1615	my $temp_std = new Statistics::Basic::StdDev(\@array);
1616	$std_dev = $temp_std -> query;
1617
1618	return $std_dev;
1619	}
1620
1621	######################################################################
1622	### Max ##############################################################
1623	sub max {
1624	#From perldoc.com
1625	my $max = shift(@_);
1626	foreach $foo (@_) {
1627	$max = $foo if $max < $foo;
1628	}
1629	return $max;
1630	}
1631
1632	######################################################################
1633	### Modulo ###########################################################
1634	sub modulo($element,$modulus) {
1635	#Do modulus of the element, but don't throw away the remainder
1636
1637	my($element,$modulus,$el,$floor) = @_;
1638
1639	if ($element < 0) {
1640	return $element + $modulus;
1641	}
1642	elsif ($element <= $modulus) {
1643	return $element;
1644	}
1645	elsif ($element > 3*$modulus) {
1646	#print "Wow $element is HUGE!\n";
1647	return $fill_value;
1648	}
1649	elsif ($element > 2*$modulus) {
1650	return $element - 2*$modulus;
1651	}
1652	elsif ($element > $modulus) {
1653	return $element - $modulus;
1654	}
1655	}
1656
1657
1658
1659

Note: See TracBrowser for help on using the browser.

root/tower_buoy/trunk/bin/tower_parse_dev.perl

Download in other formats: