NCCOOS Trac Projects: Top | Web | Platforms | Processing | Viz | Sprints | Sandbox | (Wind)

root/tower_buoy/trunk/bin/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_u**2 + $wind_v**2)**(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.