Revision Log

tower_parse.perl

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

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

Note: See TracBrowser for help on using the browser.

root/tower_buoy/trunk/bin/tower_parse.perl

Download in other formats: