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

root/tower_buoy/bin/tower_parse_dev.perl

Revision 7 (checked in by cbc, 18 years ago)

--

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