#!/usr/bin/perl use strict; use SOAP::Lite; use Time::Local; use CGI; use Chart::Graph::Gnuplot qw(gnuplot); use XML::RSS; my $verbose =1; my $query= new CGI; if ($ARGV[0] =~ /-[hv]/) { shift; $verbose++; } my $station_id = $ARGV[0] || $query->param('station') || '8637689'; # Yorktown my $startString= $ARGV[1] || $query->param('start') || 'now'; my $endString = $ARGV[2] || $query->param('end') || 48; my $scale = $query->param('scale') || 1; my $offset = $query->param('offset') || 0; my $lambda = $query->param('lambda') || .2; my $decay = $query->param('decay') || 0.997116; # 0.5 ^ (1/240) or 1 day half-life my $t0 = $query->param('t0') || ''; my $above = $query->param('above'); undef $above if defined $above && !($above=~ m/[0-9.]+/); my $below = $query->param('below'); undef $below if defined $below && !($below=~ m/[0-9.]+/); my $doRSS = $query->param('rss') || 0; my $doKML = ($query->param('kml') && ! $doRSS) || 0; my $plotfile = "tmp/coopspred$$.png"; my $startSecs = time()-3600*24 ; $startSecs -= $startSecs %3600; if ( ! ($startString =~ /^now$/)) { # +/-delta mode if ( $startString =~ m/^[-+]/){ $startSecs += $startString *3600; } else { $startString =~ m/(\d{4})-?(\d{2})-?(\d{2})[T ]?(\d{2})?:?(\d{2})?:?([0-9.]{2,})?/; my($year, $month, $day, $hour, $min, $sec )= ($1, $2, $3, $4?$4:0,$5?$5:0,$6?$6:0); # print "ISO mode $year-$month-$day $hour:$min:$sec\n"; $startSecs = timegm($sec,$min,$hour, $day, $month-1, $year-1900) ; $startSecs -= $startSecs %360; # floor 6min } } my $endSecs = $startSecs + 3600 * $endString; my $start_time = coopsTime($startSecs); my $end_time = coopsTime($endSecs); #print "$startString, $endString, $startSecs, $endSecs, $start_time, $end_time \n"; $start_time =~ s/T/ /; $end_time =~ s/T/ /; # Pass the WSDL url to the service method. But add ->want_som(1) #so we can check resonse->fault and use XPath query to loop through # response. my $soap = SOAP::Lite ->new(); my $service = $soap->service('http://opendap.co-ops.nos.noaa.gov/axis/services/Predictions?wsdl')->want_som(1); my $adjustString = ($scale != 1 || $offset != 0 )? "(Adjusted by * $scale + $offset )" :''; my $title ="Water Level $endString-hour Forecast for COOPS Station $station_id$adjustString."; $title =~ s/Forecast/Forecast above $above/ if defined $above; $title =~ s/Forecast/Forecast below $below/ if defined $below; $title =~ s/above/or above/ if ( defined $below && defined $above ); $query->delete('rss'); my $url= $query->self_url; my $rssfeed = $url . "&rss=y"; $rssfeed =~ s/t0=[^;&]*[;&]// ; # Remove t0 from the string $rssfeed =~ s/start=[^;&]*[;&]// ; # Remove start from the string #$rssfeed =~ s/end=[^;&]*[;&]/end=24;/ ; # Remove t0 from the string my $kmlfeed = $url . (($url =~ m/\?/ )? '' : '?') . '&kml=y&'; $kmlfeed =~ s/scale=[^;&?]*[;&?]?//; $kmlfeed =~ s/offset=[^;&?]*[;&?]?//; #start getting data: my $predictionParamHelp =<getPredictions($station_id, $start_time, $end_time, 0, 0,0, 6); if($response->fault()) { print 'Prediction Errors: -' . $response->faultcode . '- ' . $response->faultstring . "\n"; print $response->faultstring . " for station $station_id\n"; # exit; } # Pass the WSDL url to the service method. But add ->want_som(1) so we can check resonse->fault and use XPath query to loop thru # response. my $service = SOAP::Lite->service('http://opendap.co-ops.nos.noaa.gov/axis/services/WaterLevelRawSixMin?wsdl')->want_som(1); my $obs = $service->getWaterLevelRawSixMin($station_id, $start_time, $end_time, 'MLLW', 0, 0); if($obs->fault()) { print 'Error: -' . $obs->faultcode . '- ' . $obs->faultstring . "\n"; print $obs->faultstring . " for station $station_id\n"; # exit; } #print "$station_id\n"; #for my $t ($obs->valueof('//data/item') ) { # print "$t->{timeStamp} $t->{WL}\n"; # # some QA flags. # #print "$t->{sigma} $t->{O} $t->{F} $t->{R} $t->{L}\n"; #} # #print "$station_id\n"; #for my $t ($response->valueof('//data/item') ) { # # # COOPS prediction timestamps are formatted different than observations: # my ($mon, $day, $year, $time) = split '/| ',$t->{timeStamp}; ## my $hashlabel = $t->{timeStamp} . ':00.0'; # my $hashlabel = sprintf("%04d-%02d-%02d $time:00.0",$year,$mon,$day,$time); # print "$hashlabel $t->{pred}\n"; #} my @preds = ($response->valueof('//data/item')); my @obss = ($obs->valueof('//data/item') ); my $datastring ="
\n";

$datastring .= "timestamp(UTC)    predict   obs  resid  sm_resid forecast, scaled\n" if $verbose;
my $smresidual = 0;
my $residual=0;
my @predstamp;
my @predvals;
my @obsvals;
my @resids;
my @smresids;
my @forecasts;
my $maxforecastii = '';
my $minforecastii = '';
my $lastobsii;
my $ij=0 ;
for (my $ii = 0 ; $ii <= $#preds ; $ii++) {
#    print "$ii, ", $preds[$ii]->{timeStamp};
#    print " ", ($obss[$ii]->{timeStamp} || " . "), ;
#    print $preds[$ii]->{pred}, " ", ($obss[$ii]->{WL} || " . ");
#    print "\n";
    my ($month,$day,$year,$time) = split '[/ ]',$preds[$ii]->{timeStamp};
    push @predstamp,sprintf("%04d-%02d-%02dT%s",
			    $year,$month,$day,$time);
    my $pred = $preds[$ii]->{pred};
    my $x = -999;
    $residual ='x';
    my $forecast = $pred + $smresidual;
    my $scForecast = $forecast *$scale + $offset;
    $maxforecastii = $ii if ( ($maxforecastii eq '' ) || 
			      $scForecast > $forecasts[$maxforecastii]);
    $minforecastii = $ii if ( ($minforecastii eq '') || 
			      $scForecast < $forecasts[$minforecastii]);
# bad assumption-- make it so missing data is aligned
#    if (  $obss[$ii]->{WL} ) {  # use new data for updating error
    if (  $obss[$ij]->{timeStamp} eq sprintf('%04d-%02d-%02d %s:00.0',$year,$month,$day,$time)) { 
 # use new data for updating error
	$x = $obss[$ij]->{WL};
        $residual = $x-$pred;
        if ( ! $t0 || $t0 ge $predstamp[$ii]){
	    $smresidual = $smresidual ? 
		$smresidual *(1-$lambda) + $lambda*($residual ) 
		: $residual;
	}
	$ij++;
	$lastobsii = $ii;
    }

    $datastring .= sprintf("%04d-%02d-%02dT%s %7.3f %8.3f %7.3f %7.3f %7.3f %7.3f\n",
          $year,$month,$day,$time, # $obss[$ij]->{timeStamp} || '9999-99-99 99:99:99',
    $pred , 
    $x, 
    $residual,
    $smresidual, 
    $forecast, $forecast *$scale +$offset);

    $smresidual=$smresidual * $decay ; #0.997116; # 0.5 ^ (1/240) or 1 day half-life    

    push @predvals, $pred * $scale + $offset;
    push @resids, $residual * $scale;
    push @smresids, $smresidual * $scale;
    push @obsvals, $x == -999 ? 'x':$x * $scale + $offset ;
    push @forecasts, $scForecast;

}

$datastring .= "
\n"; my $minWL = sprintf("%.3f",$forecasts[$minforecastii]); my $maxWL = sprintf("%.3f",$forecasts[$maxforecastii]); my $nowWL = sprintf("%.3f",$forecasts[$lastobsii]); $| = 1; if ( $doRSS ) { doRSS(); } elsif ( $doKML ) { doKML(); } else { print $query->header; print $query->start_html(-title =>$title, -head =>[$query->Link({-rel=>'alternate', -type=>"application/rss+xml", -title=>"$title RSS", -href=>"$rssfeed", }), $query->Link({-rel=>'alternate', -type=>"application/vnd.google-earth.kml+xml", -title=>"$title KML", -href=>"$kmlfeed", }), ], -meta => {'robots' => 'nofollow'}, ); doPlot(); sleep 2; print qq!

$endString-hour EWMA forecast of Station $station_id $adjustString

Max/Min water level forecasts: $maxWL / $minWL
!; print_prompt($query); print_tail(); print '

Data

Data derived from COOPS SOAP services:WSDL(observations) WSDL(tides)',$datastring; print $query->end_html; #open( STDOUT,">/dev/null"); #open( STDERR,">/dev/null"); #sleep 60; #unlink $plotfile; } BEGIN { sub coopsTime { my $utcSecs = shift; my ($sec,$min,$hour,$day,$mon0,$yr1900,$dow,$isdst) = gmtime($utcSecs); # wrfout_d01_2006-08-30_00:00:00.out my $format= "%04d%02d%02d %02d:%02d"; return sprintf($format,$yr1900+1900, $mon0+1, $day, $hour, $min); } sub doKML{ #Content-Type: application/vnd.google-earth.kml+xml #Content-Type: text/xml print < ${maxWL}m Max Station $station_id flood surface 1 0 This is a simple plane $maxWL m above MSL at its corners. See $url for more info. 0 0 http://sura-vims-pe6600-1.vims.edu/~drf/SCOOP/viewCenteredFlood.cgi?FLOOD=$maxWL 2 onRequest 1 EOT } sub doKMLi{ #Content-Type: application/vnd.google-earth.kml+xml #Content-Type: text/xml print < ${maxWL}m flood surface -122.0822035425683,37.42228990140251,0 EOT } sub doRSS{ my $rss = new XML::RSS (version => '2.0'); $rss->channel(title => $title, link => "$url", language => 'en', description => "Surge-corrected Tidal predictions for COOPS station $station_id $adjustString", copyright => 'Copyright 2006, David Forrest', pubDate => 'Mon Nov 20 16:05:11 2006', lastBuildDate => sprintf("%s",my $dummy=gmtime()), docs => 'http://sura-vims-pe6600-1.vims.edu/~drf/coopspred.pl', managingEditor => 'drf@vims.edu', webMaster => 'drf@vims.edu' ); $rss->add_module(prefix=>'georss', uri=>'http://www.georss.org/georss'); $rss->add_module(prefix=>'gml', uri=>'http://www.opengis.net/gml'); if (( (! defined $above && ! defined $below ) || # no threshholding ( defined $above && $maxWL > $above) || # too high ( defined $below && $minWL < $below) )) { #too low $rss->add_item(title => sprintf("COOPS %s: %7.3f [%7.3f] %7.3f / ${endString}hr", $station_id, $minWL, $forecasts[$lastobsii], $maxWL, $endString), # creates a guid field with permaLink=true # permaLink => "http://freshmeat.net/news/1999/06/21/930003829.html", # alternately creates a guid field with permaLink=false guid => "$url;rsst0=$predstamp[$lastobsii]", link => "$url", description => <<"EOT", $title The COOPS station $station_id reading $adjustString is now $nowWL and may vary between $minWL and $maxWL over the next $endString hours. EOT georss => {point =>"0 0"}, # GeoRSS simple works # next lines fail by not expanding the gml hashes: # i georss => {point =>"0 0", # where => gml => { Point => # gml => { pos => "0 0"}}, # }, # Next line doesn't work because it escapes the XML: # georss => { where => "0 0",} ) } else { my $threshhold = (defined $above) ?"$above high limit":''; $threshhold =~ s/$/ or / if (defined $above && defined $below); $threshhold =~ s/$/$below lower limit/ if (defined $below); $rss->add_item(title => sprintf("OK: COOPS %s: %7.3f [%7.3f] %7.3f / ${endString}hr", $station_id, $minWL, $forecasts[$lastobsii], $maxWL, $endString), # creates a guid field with permaLink=true # permaLink => "http://freshmeat.net/news/1999/06/21/930003829.html", # alternately creates a guid field with permaLink=false guid => "$url;", link => "$url", description => <<"EOT", $title The COOPS station $station_id reading $adjustString is now $nowWL and may vary between $minWL and $maxWL over the next $endString hours. However, this model does not expect violation the threshhold of $threshhold. EOT ) } my $data= "Content-Type: text/xml\n\n" . $rss->as_string; my $stylesheet=''; my $xslsheet=''; $data =~ s/(<\?xml[^>]+>)/$1$stylesheet$xslsheet/; print $data; } sub doPlot { mkdir 'tmp',0755 if ( ! -d 'tmp' ) ; gnuplot({'title' => "Water level for Station $station_id $adjustString", 'x-axis label' => 'Date and Time(UTC)', 'y-axis label' => "$scale X (Meters above $station_id MLLW)+ $offset", 'output type' => 'png', 'output file' => $plotfile, # Setting date/time specific options. 'xdata' => 'time', 'timefmt' => '%Y-%m-%dT%H:%M', 'format' => ['x', '%Y/%m/%d\n%H:%M'], # Set output range - note quoting of date string # 'xrange' => '["' . $predstamp[1] .'":"' . $predstamp[50].'"]', 'extra_opts' => join("\n", 'set grid', 'set timestamp'), }, # Data [{'title' => "Observation at $station_id", 'type' => 'columns', 'style' => 'points', }, [ @predstamp ], [@obsvals], ], # Data [{'title' => 'Astronomical Tidal Prediction', 'type' => 'columns', 'style' => 'lines', }, [ @predstamp], [@predvals], ], # Data [{'title' => 'Residual', 'type' => 'columns', 'style' => 'lines', }, [@predstamp], [@resids], ], # Data [{'title' => 'EWMA Smoothed Residual', 'type' => 'columns', 'style' => 'lines', }, [@predstamp], [@smresids], ], [{'title' => 'Surge-based EWMA Forecast', 'type' => 'columns', 'style' => 'lines', 'axes' => 'x1y1', }, [@predstamp], [@forecasts], ], [{'title' => "Max. Forecast: $maxWL", 'type' => 'function', 'style' => 'lines', }, "$maxWL", ], [{'title' => "Min. Forecast: $minWL", 'type' => 'function', 'style' => 'lines', }, "$minWL", ], # ( (defined $above) ? [{'title' => "Max. Threshhold: $above", # 'type' => 'function', # 'style' => 'lines', # }, # "$above", # ]:''), # # ( (defined $below) ? [{'title' => "Min. Threshhold: $below", # 'type' => 'function', # 'style' => 'lines', # }, # "$below", # ]:''), ); } sub print_prompt { my($query) = @_; print $query->start_form(-method=>'GET'); print 'COOPS station?'; print $query->textfield('station',$station_id,'size',7); print "Start:"; print $query->textfield('start',$startString); print "End:"; print $query->textfield('end',$endString); print "
\nscale:"; print $query->textfield('scale',$query->param('scale')||1); print "offset:"; print $query->textfield('offset',$query->param('offset')||0); print "
\nAdvanced Options:
\n"; print "t0 (for retrospective forecasts):", $query->textfield('t0',$query->param('t0')|''); print "
\nlambda (for different EWMA smoothing):", $query->textfield('lambda',$query->param('lambda')|'0.2'); print "
\ndecay (for different EWMA forecast):", $query->textfield('decay',$query->param('decay')|'0.997116'), " 1.0 for clamp, 0.997116 for 1-day halflife"; print "
Threshhold above:", $query->textfield('above',$query->param('above')|''), qq! (upper limit for RSS feed)!; print "
Threshhold below:", $query->textfield('below',$query->param('below')|''), qq! (lower limit for RSS feed)!; print "

",$query->reset; print $query->submit('Action','Submit'); print $query->endform; print "

\n"; } sub do_settings { my($query) = @_; my(@values,$key); print "

Here are the current settings in this form

"; foreach $key ($query->param) { print "$key -> "; @values = $query->param($key); print join(", ",@values),"
\n"; } } sub print_tail { print "Max(forecast)= ",sprintf("%7.3f at %s
\n", $maxWL, ,$predstamp[$maxforecastii]); print "Min(forecast)= ",sprintf("%7.3f at %s
\n", $minWL, ,$predstamp[$minforecastii]); print <Data from $station_id at COOPS, Historic data in meters. For more information on this station, see NOAA Tides And Currents Station $station_id. This tool uses the NOAA COOPS Web Services, as can be seen in the source code.

How to use this model:

  • Choose a COOPS station close to your area of interest and plug it into the 'station' box above. By default, the program will train itself on the last 6 hours of observations, and use that to estimate surge for the following 24 hours.
  • You can try to calibrate the station to a more local area of interest (the end of your dock, for instance) by scaling and shifting the COOPS observation to match. As an example, Greate road near VIMS is about 5\' above the Yorktown station 8637689, so scale the meters by 3.2808 and add 5\' to get water over the Greate Road.
  • You can look at historical data by choosing a different start time, either as an offset from now with a -hhhhh where hhhhh is the number of hours to go back in time, or with an absolute start time in the form yyyy-mm-ddThh:mm.
  • The t0 setting limits the training set to end at t0, so that you can see what would have been known and predicted by this model at a certain time in the past. 't0' should also be in the form 'yyyy-mm-ddThh:mm'.
  • There are a couple pre-configured examples below.
  • The raw and derived data is included below.

Model Details

Forecast model: Forecast = tides + forecast residuals, where the forecast residuals are calculated using EWMA smoothing: lambda=.2 retrospective, 1-(0.5 ^ (1/240))=0.00288 exponential decay of the forecast residuals.

The EWMA smoothing on input is intended to smooth out some of the noise in the observation -- the smoothed estimate of the residual is 20% due to the most recent observation, and 80% due to the previous estimate. If you would like to try a different smoothing coefficient, try using the parameter 'lambda=0.8'.

Bugs: The code writes its graphics to a temporary filename ($plotfile), which may disappear by the time your browser downloads this page.

[RSS] newsfeed for $title. Add it as a live bookmark in Firefox, or a MS IE7 feed.

[KML] For the maximum water level over as a polygon/building in googleEarth -- Take care with this as the datums are not consistent.Use offset=(msl-mllw) to get to MSL, which is what GE says it uses. David Forrest 2007-03-13

The prediction model used here is a simple assumption that the surge in the future will be similar, but lower than what we see right now. It does not take into account any effects of weather water, precipitation, or stream flow dynamics. VIMS is working on a hydrodynamic model for surge forecasting which takes some of these variables into account.

Inputs & Outputs:

  • Station: A COOPS station number that works with ther SOAP services
  • Start: An offset in hours for the beginning time of the data collection. It also works with an absolute ISO-8601 timestamp like 2005-08-02T18:00 or the special word 'now' to start taking data from 1 hour ago.
  • End: An offset in hours for the end time of the data collection. It also works with an absolute ISO-8601 timestamp like 2005-08-02T18:00.
  • scale: a multiplier to change the forecast from meters to any other unit
  • offset: a constant to add to each forecast to allow for calibration to other places and gauges
  • Timestamp: the UTC time of the particular data point. UTC is 4 hours ahead of EDT.
  • predict: The COOPS predction from tidal constituents at the station (meters)
  • obs: The observed water level at the station (m)
  • sm_resid: Am EWMA smoothed residual of the observation (m). In forecast mode, the forecast residual decays exponentially with a half-life of 1 day.
  • forecast: The COOPS tidal prediction+the prior smoothed residual (m)
  • scaled: the forecast*scale+offset for calibrating forecasts to other locations.

Examples:

EOT print <David Forrest -- 2006-10-06
END } }