#!/afs/isis/pkg/isis/bin/perl
# Last modified:  Time-stamp: <2003-11-06 12:15:01 haines>
#
# Abstract:  get cycle metar data from NWS, and archive the file based on 
#    ICAO station ids to raw data directory
#
# Usage:  % archive_cycle_metar.pl [-d]
#
# Author: Sara Haines (2003-06-19) 
#

# processing steps
# (1) determine which METAR cycle file to get based on current time (UTC)
# (2) read list of stations to determine which stations 
# (3) for each station, 
#      a. find latest obs in cycle file
#      b. determine if latest obs is unique to archive data
#      c. write data to output file 

$delimiter = "\n";

# if debugging requested print messages to STDOUT
if (grep /[debug|DEBUG|d]/, @ARGV) {
  $debug = 1;
}

# add libraries needed for this function
use POSIX qw(strftime);
use POSIX qw(strtod);
use POSIX qw(exp);
use POSIX qw(sqrt);
use CGI qw(:standard);

use Date::Manip;
use Time::Local;
use LWP::Simple;

# print Date and Time and version of perl + patchlevel /1000
# way to do date/time manip through Date::Manip
$now=&UnixDate(&ParseDate("today"), "%Y:%m:%d %H:%M:%S");
if ($debug) { print "\n==== Starting: $now ==== Perl Version: $]\n"; }

# raw data paths
$raw_path = "/seacoos/data/nws_metar/raw_data";

# some date things
$this_month=&UnixDate(&ParseDate("today"), "%Y_%m");
$prev_month=&UnixDate(&DateCalc("today", "- 1 month"), "%Y_%m");
$this_day = (gmtime)[3];

# start stop-watch for cummulative times
$starttime = time;


# -----------------------------------------------------------------------
# (2) read list of stations to determine which stations
# -----------------------------------------------------------------------
# list of station ids to search for in metar cycle file
# just create it for now
# @station_list = ("KRDU", "KIGX", "KTTA", "KSAV", "KILM");
# @station_list = ("KRDU");

$doc = get("file:/opt/local/seacoos/bin/stations_seacoos.txt");

# splitting a new line drops new lines.
@station_data = split(/$delimiter/, $doc);

# if (debug) { print $doc;}

# extract station_ids
# there's probably a better, cleaner way to do this in perl, but ...
foreach $_ (@station_data) {
	push @station_list, substr($_, 0, 4);    
}

# ****
# @station_list = ("KRDU", "KIGX", "KTTA", "KSAV", "KILM");

# -------------------------------------------------------$station_dir ----------------
# (3) for each station, 
#      a. find latest obs in cycle file
#      b. determine if latest obs is unique to archive data
#      c. write data to current month output file 
# -----------------------------------------------------------------------

$starttime2 = time;

for ($i=0; $i<=$#station_list; $i++) {

    $line1 = "";
    $station_id = $station_list[$i];
    $pattern = "(" . $station_list[$i] . ")";
    
    # get station METAR
    $station_file = "http://weather.noaa.gov/pub/data/observations/metar/stations/" . 
	$station_id . ".TXT";
    $doc = get($station_file);

    # splitting a new line drops new lines.
    @metar_data = split(/$delimiter/, $doc); 

    # find lines in @metar_data that contains station id
    # should only be one line, but just in case there is more than one
    @metar_data_subset = grep /$pattern/, @metar_data;

    if ($#metar_data_subset == -1) {
	$have_data = 0;
    } else {
	# get last metar found 
	# last one should be last recorded obs for that cycle
	$line1 = $metar_data_subset[$#metar_data_subset];
	if ($line1) { $have_data = 1; }
    }

    if ($have_data) {

	# determine if station has a archive directory yet
	$station_dir = ($raw_path . "/" . $station_id);

	# if station dir does not exist, then mkdir 
	if ( !(-e "$station_dir") ) {
	    `mkdir $station_dir`;
	}

	# Handle month transition and latency of data
	# e.g. current time is 2003-09-01 01:00 but data is for hour 23 (hence 2003-08-31)

	# parse day of month from coded metar ($DD)
	&parse_metar_time; 

	# compare to current day of month ($DD_current) to determine which 
	# archive file to use, previous month or current month
	# current month (cm) file
	if ($this_day >= $DD) {
	    $cmfile = ($station_dir . "/" . $station_id . "_met_" . $this_month . ".dat");

	} else {

	    $cmfile = ($station_dir . "/" . $station_id . "_met_" . $prev_month . ".dat");
	}
	
	# if (debug) { print($cmfile, "\n"); }

	if ( !(-e "$cmfile") ) {
	    # if data file does not exist
	    # output $line1 to new file
	    open(FILE, ">$cmfile");
	    print FILE "$line1\n";
	    close(FILE);

	    # if (debug) { print($line1, "\n"); }
	    
	} else {
	    # if data file does exists

	    # open and read entire file input and stores it as an array
	    open(FILE, "$cmfile");
	    @FILE = <FILE>;
	    close(FILE);

	    # joining with a null and splitting a new line effectively drops new lines.
	    $chunks = join('',@FILE);
	    @month_data = split(/$delimiter/, $chunks);
	    
	    # determine if $line1 matches any line in @month_data
	    # first replace any regexp char with escaped char
	    $line1_str = $line1;
	    $line1_str =~ s/([\*\.\<\>\{\}\[\]\^\$\|\+\?\\\/])/\\$1/g;
	    # print "$line1\n";
	    $line_pattern = "(" . $line1_str . ")";
	    @match = grep /$line_pattern/, @month_data;

	    if ($#match == -1) {
		# if there are no matches, can assume its a 
		# new observation and need to append to current month data

		open(FILE, ">>$cmfile");
		print FILE "$line1\n";
		close(FILE);

		# if (debug) { print($line1, "\n"); }
	    }
	}
    }


    # &parse_metar;
    # &print_full_text;
    # &print_line_text;
    # &reset_obs;

} # for ($i=0; $i<=$#station_list; $i++)

$subtime = time - $starttime2;
if ($debug) { print (" ... archive time = $subtime (seconds)\n"); }

$cummtime = time - $starttime;
if ($debug) { print ("Total script time = $cummtime (seconds)\n"); }
# if ($debug) { print("... Done\n"); }

###############################subroutines ##################################
sub parse_metar_time {

# 
# Major parts and general idea of this subroutine courtesy of: 
# wmk <wkowalsk@ii.pw.edu.pl>
# www.ii.pw.edu.pl/~wkowalsk/metar2.cgi
# Thu May 20 21:34:39 CET 1999
#

    $data = $line1;
    $data =~ s/^($station_id )//g;
    @pos = split(/ /,$data);

    $rtime = "";

    foreach $_ (@pos) {
	if    (/RMK|TEMPO|BECMG|NOSIG/){
	    goto NEXT; 
	}
	elsif (/^(\d{4})Z|^(\d{6})Z|^(\d{6})/){	
	    if (!$rtime){
		$rtime = $_;	
		goto NEXT; 
	    }
	}
    } #  foreach $_ (@pos)

  NEXT:

    @times = split (//,$rtime);
    if (@times > 5)
    {
	($t1,$t2,$t3,$t4,$t5,$t6) = @times;
    }
    else
    {
	($t3,$t4,$t5,$t6) = @times;
    }

  END:

    if (!$rtime) {
	$t1 = strftime("%d",gmtime);
	$t2 = ''; 	
	# $t1 = strftime("%d %H %M",gmtime);
	# $t3 = '';$t2 = '';$t4 = '';$t5 = ''; $t6 = '';
    }
    elsif (!$t1 && !$t2) {
	$t1 = strftime("%d",gmtime);
	$t2 = ''; }

    # if (debug) { print "as of $date $t1$t2 $t3$t4 $t5$t6 UT\n"; }
    $DD = join("", $t1, $t2);
}


sub parse_metar {

# 
# Major parts and general idea of this subroutine courtesy of: 
# wmk <wkowalsk@ii.pw.edu.pl>
# www.ii.pw.edu.pl/~wkowalsk/metar2.cgi
# Thu May 20 21:34:39 CET 1999
#

    $data = $line1;
    $line1 =~ s/^($station_id )//g;
    @pos = split(/ /,$line1);


    foreach $_ (@pos) {
	if    (/RMK|TEMPO|BECMG|NOSIG/){
	    goto NEXT; 
	}
	elsif (/^(\d{4})Z|^(\d{6})Z|^(\d{6})/){	
	    if (!$rtime){
		$rtime = $_;	
	    }
	}
	elsif (/(..)MPS|(..)KT/){
	    $wind1 = $_;
	}
	elsif (/^(\d{4})|(.)SM|CAVOK/){			
	    if (!$visib){
		$visib = $_;    
	    }
	}
	elsif (/BR|DZ|DS|DU|RA|SN|SG|SA|SS|SQ|IC|PE|PY|GR|GS|FG|FC|FU|HZ/){	
	    push @phenom, $_;
	}
	elsif (/FEW|SCT|BKN|OVC/){	
	    push @clouds, $_;
	}
	elsif (/^(\d{3})V(\d{3})/){
	    if (!$wind2){
		$wind2 = $_; 
	    }
	}
	elsif (/^(M?\d{2})\/(M?\d{2})/){
	    if (!$tempr){	
		$tempr = $_;
	    }
	}
	elsif (/^Q(\d{4})|^A(\d{4})/){
	    if (!$press) {
		$press = $_;
		goto NEXT;
	    }
	}
    }


 NEXT:

#    # search for stuff beyond RMK
#     foreach $_ ($pos) {
# 	if (/^T/) {
# 	    # $tempr =~ s/M/\-/g;
# 	    # $tempr = $_;
# 	    # $temp_air = substr($tempr, 1,4)/10;
#  	    # $temp_dew = substr($tempr, 5,4)/10;
# 	    print "HAVE HIGHER PRECISION temperature\n";
# 	    # print ($temp_air, " ", $temp_dew, "\n");
# 	}
# 	if (/^SLP/) {
# 	    print "HAVE HIGHER PRECISION pressure\n";
# 	}
#     }

    $visib =~ s/CAVOK|9999|10SM/Unlimited/g;
    if ($visib =~ s/SM//g) {
	$visib *= 1600; 
    }
    if ($visib ne 'Unlimited') {
	$visib *= 0.001;
	$visib = join("", $visib, " km"); 
    }


    @times = split (//,$rtime);
    if (@times > 5) {
	($t1,$t2,$t3,$t4,$t5,$t6) = @times;
    }
    else {
	($t3,$t4,$t5,$t6) = @times;
    }


    $press =~ s/Q//g;
    if ($press =~ s/A//g) {
	$press *= 0.33868;
	$press = sprintf("%.f", $press);
    }

    $tempr =~ s/M/\-/g;
	($temp_air,$temp_dew) = split(/\//,$tempr);

    if ($wind2 =~ s/V/-/g) {
	$wind2 = join("", ", variability: ", $wind2, " deg");
    }
    if ($wind1 =~ m/00000/) {
	$speed = "Calm";
    }
    else {
	$mps = "m/s";
	$wind1 =~ s/MPS//g;
	if ($wind1 =~ s/KT//g){	$kt = 1; }
	if ($wind1 =~ m/VRB/){	$vr = 1; }
	@winds = split (//,$wind1);
	if (@winds > 7) {
	    ($w1,$w2,$w3,$w4,$w5,$w9,$w10,$w11) = @winds;
	    $w9 =~ s/G/, gusting to /g;
	    $speed = join("", $w4, $w5);
	    $speeg = join("", $w10, $w11);
	    if ($kt) {
		$speed = sprintf("%.f", $speed*0.514);
		$speeg = sprintf("%.f", $speeg*0.514);
	    }
	    # $speeg = strtod($speeg);
	    $speeg = join(" ", $speeg, $mps);
	} 	
	else {	
	    ($w1,$w2,$w3,$w4,$w5) = @winds;
	    $speed = join("", $w4, $w5);
	    if ($kt) {
		$speed = sprintf("%.f", $speed*0.514);
	    }
	}
	if (!$vr) {
	    $wdir = join ("", $w1,$w2,$w3);
	    $w1 = join("", "from ", $w1);
	    $w3 = join("", $w3, " deg");
	    
	    if    ($wdir < 15){	$wdir_str = ' (N)';  }
	    elsif ($wdir < 30){	$wdir_str = ' (NNE)';}
	    elsif ($wdir < 60){	$wdir_str = ' (NE)'; }
	    elsif ($wdir < 75){	$wdir_str = ' (ENE)';}
	    elsif ($wdir < 105){$wdir_str = ' (E)';  }
	    elsif ($wdir < 120){$wdir_str = ' (ESE)';}
	    elsif ($wdir < 150){$wdir_str = ' (SE)'; }
	    elsif ($wdir < 165){$wdir_str = ' (SSE)';}
	    elsif ($wdir < 195){$wdir_str = ' (S)';  }
	    elsif ($wdir < 210){$wdir_str = ' (SSE)';}
	    elsif ($wdir < 240){$wdir_str = ' (SW)'; }
	    elsif ($wdir < 265){$wdir_str = ' (SSW)';}
	    elsif ($wdir < 285){$wdir_str = ' (W)';  }
	    elsif ($wdir < 300){$wdir_str = ' (WNW)';}
	    elsif ($wdir < 330){$wdir_str = ' (NW)'; }
	    elsif ($wdir < 345){$wdir_str = ' (NNW)';}
	    else  { $wdir_str = ' (N)';	}
	}
	else	{
	$w1 = "direction variable";
	$w2 = '';$w3 = '';
    }


	$speew = $speed;
	# $speed = strtod($speed);

	if ($speew > 2)
	{
	    $refroid = ((((sqrt($speew*100))+10.45-$speew)*(33-$temp_air))*1.163);
	    $temp_wch = sprintf("%.f", (33-$refroid/26.905871));
	}
    }

    $esdb = exp((($temp_air*16.78)-116.9)/($temp_air+237.3));
    $eswb = exp((($temp_dew*16.78)-116.9)/($temp_dew+237.3));
    $humid = sprintf("%.f",(($eswb/$esdb)*100));


    foreach $_ (@clouds) {
	~ s/([0-9]+)/($1*0.030)/eg;
	~ s/FEW/Few /g;
	~ s/SCT/Scattered /g;
	~ s/BKN/Broken /g;
	~ s/OVC/Overcast /g;
	~ s/TCU/ Towering Cumulus/g;
	~ s/TCB/ Towering Cumulonimbus/g;
	~ s/CU/ Cumulus/g;
	~ s/CB/ Cumulonimbus/g;
    }
    foreach $_ (@phenom) {
	~ s/DR/Drifting /g;
	~ s/MI/Shallow /g;
	~ s/PR/Partial /g;
	~ s/FZ/Freezing /g;
	~ s/BL/Blowing /g;
	~ s/SH/Shower /g;
	~ s/BC/Patches /g;
	~ s/TS/Thunderstorm /g;
	~ s/\+FC/Tornado/g;
	~ s/FC/Funnel Cloud/g;
	~ s/DS/Duststorm/g;
	~ s/DU/Widespread Dust/g;
	~ s/SA/Sand/g;
	~ s/SS/Sandstorm/g;
	~ s/SQ/Squall/g;
	~ s/RASN/Rain and Snow/g;
	~ s/SNRA/Rain and Snow/g;
	~ s/SN/Snow/g;
	~ s/PY/Spray/g;
	~ s/FU/Smoke/g;
	~ s/GR/Hail/g;
	~ s/GS/Small Hail/g;
	~ s/RA/Rain/g;
	~ s/DZ/Drizzle/g;
	~ s/BR/Mist/g;
	~ s/SG/Snow Grains/g;
	~ s/IC/Ice Crystals/g;
	~ s/PE/Ice Pellets/g;
	~ s/FG/Fog/g;
	~ s/HZ/Haze/g;
	~ s/UP/*Unknown*/g;
	~ s/\-/Light /g;
	~ s/\+/Heavy /g; }
    if (@phenom > 1) {
	@phenom = join(", ", @phenom); }
    if (@clouds >= 1) {
	@clouds = join(", ", @clouds);
	push @clouds, '[km]'; }


 END:

    if (!$rtime) {
	$t1 = strftime("%d %H %M",gmtime);
	$t3 = '';$t2 = '';$t4 = '';$t5 = ''; $t6 = ''; }
    elsif (!$t1 && !$t2) {
	$t1 = strftime("%d",gmtime);
	$t2 = ''; }

    $date = strftime("%Y %m",gmtime);

    #&print_full_text;
}

# ===================== output different formats ====================

sub print_full_text {

    print "-----------------------------------------\n";
    print "Weather Conditions at $station_id ";
    print "as of $date $t1$t2 $t3$t4 $t5$t6 UT:\n";

    if ($nodata) {
	print "No data available\n";
    }

    else {
	print "Wind: $speed $mps$w9$speeg $w1$w2$w3$wdir_str$wind2\n";
	print "Visibility: $visib\n";
	if (@clouds){	print "Clouds: @clouds\n"; }
	print "Temperature: " . strtod($temp_air) . " deg C\n";
	if ($temp_wch){	print "Windchill: " . strtod($temp_wch) . " deg C\n"; }
	print "Dew point: " . strtod($temp_dew) . " deg C\n";
	print "Relative Humidity: $humid %\n";
	print "Pressure: " . strtod($press) . " hPa\n";
	if (@phenom){	print "Weather: @phenom\n"; }
    }


}

sub print_line_text {

    if ($nodata) {
	# do nothing
    }
    else {
	# print ICAO STATION CD LAT LON ELEV 
	print "$station_data[$i] ";
	# print YYYY MM DD hh mm 
	print "$date $t1$t2 $t3$t4 $t5$t6 ";
	# WDIR WSPD
	if ($vr) { print "999 $speed "; }
	elsif ($speed == 'Calm') { print "000 0 "; }
	else { print "$wdir $speed "; }
	
	# ATEMP DEWP PRESS 
	if (!$temp_air) {print "999 ";} else {print("$temp_air ");}
	if (!$temp_dew) {print "999 ";} else {print("$temp_dew ");}
	if (!$press) {print "999 ";} else {print("$press ");}
	# EOL
	print "\n";
    }
}

sub reset_obs {

    $nodata = "";
    $rtime = "";

    @winds = "";
    $wind1 = $wind2 = "";

    $vr = $kt = "";
    $speed = $mps = $w9 = $speeg = ""; 
    $w1 = $w2 = $w3 = "";
    $wdir = "";
    $wdir_str = "";


    $visib = "";
    $tempr = "";
    $temp_wch = "";
    $temp_air = "";
    $temp_dew = "";

    $press = "";

    @clouds = "";
    @phenom = "";
}