#!/usr/bin/perl ## import CGI and its HTML methods use CGI qw/:standard start_ul end_ul/; ## import Database Interface module use DBI; use DateTime; ## this is called if a fatal exception is thrown sub last_breath { print p, "Sorry, the database is currently down! ", "Please try again later."; print p, "But you may enjoy the error message if you want to: ", $_[0]; print end_html; } # # http://www.zipcodeworld.com/samples/distance.pl.html # Simple distance from lat lon calculator # #::: Hexa Software Development Center © All Rights Reserved 2004 ::: $pi = atan2(1,1) * 4; sub distance { my ($lat1, $lon1, $lat2, $lon2, $unit) = @_; my $theta = $lon1 - $lon2; my $dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta)); $dist = acos($dist); $dist = rad2deg($dist); $dist = $dist * 60 * 1.1515; if ($unit eq "K") { $dist = $dist * 1.609344; } elsif ($unit eq "N") { $dist = $dist * 0.8684; } return ($dist); } #::: This function get the arccos function using arctan function ::: sub acos { my ($rad) = @_; my $ret = atan2(sqrt(1 - $rad**2), $rad); return $ret; } #::: This function converts decimal degrees to radians ::: sub deg2rad { my ($deg) = @_; return ($deg * $pi / 180); } #::: This function converts radians to decimal degrees ::: sub rad2deg { my ($rad) = @_; return ($rad * 180 / $pi); } ## take a last breath before dying $SIG{__DIE__} = \&last_breath; ## the HTTP header print header; ## the HTML form print start_html('CFSR Query'), h1('CFSR Query'), start_form(-method=>"GET"), p,"Lat? ",textfield('lat'), p,"Lon? ",textfield('lon'), p,"Time Offset? ",textfield('timeoff',0),"used to adjust for morning or afternoon base measurements.", p,"email? ",textfield('emailaddr'), p,"Interpolation Power ",textfield('interppow',2.0), p,submit, reset, end_form, hr; ## param returns the list of supplied parameters; the empty list ## is false in a boolean context $remoteHost = $ENV{'REMOTE_HOST'}; $remoteIP = $ENV{'REMOTE_ADDR'}; print p,"Hello ",$remoteHost," at:",$remoteIP,"\n"; if(!(param('lat'))){ print end_html; exit; } $lat= param('lat'); $lon= param('lon'); $emailaddr= param('emailaddr'); $timeoff= param('timeoff'); $tzoffset=$lon/15-($timeoff); $tzoffset=sprintf("%.0f",$tzoffset); $interppow= param('interppow'); $interppow=sprintf("%.0f",$interppow); print p,"Assuming total time zone offset (Timezone and measure) for UTC:",$tzoffset,"\n"; $database="drfuka_cfsr"; $hostname="localhost"; $port=3306; $user="drfuka_cfsr"; $password="1rhQW!"; $dsn = "DBI:mysql:database=$database;host=$hostname;port=$port"; $dbh = DBI->connect($dsn, $user, $password) or die $DBI::errstr; my $start = DateTime->now; #insert into requests (lat,lon,datetime,email,offset) values (20,-20,now(),"dan@dan.com",3) $stmt="insert into requests (lat,lon,datetime,email,offset) values ($lat,$lon,now(),\"$emailaddr\",$tzoffset)"; $sth = $dbh->prepare($stmt); $sth->execute(); # # Create a MEMORY table to hold the extracted vectors aligned in the previous step. This table will be used # to calculate the HR and wind as well as to order by time the data series. # $tmptablename=$$; $stmt="create table t${tmptablename} (datetime datetime NOT NULL, lon float NOT NULL, lat float NOT NULL,tmp2m float not null, prate float NOT NULL, q2m float not null, pressfc float not null, uwnd10m float not null, vwnd10m float not null, dswsfc float not null, wind float not null, rh float not null) ENGINE=MEMORY"; #$stmt="create table t${tmptablename} (datetime datetime NOT NULL, lon float NOT NULL, lat float NOT NULL,tmp2m float not null, prate float NOT NULL, q2m float not null, pressfc float not null, uwnd10m float not null, vwnd10m float not null, dswsfc float not null, wind float not null, rh float not null) "; $sth = $dbh->prepare($stmt); #print p,"Query:",$stmt,"\n"; $sth->execute(); # # Loop over cfsr 5 year merge tables to increase performance # @cfsr_mrgs=qw(85); @cfsr_mrgs=qw(80 85 90 95 00 05); foreach $yr (@cfsr_mrgs){ # # SQL statement to extract the month of hourly vectors for each of the variables from the compressed blobs # for the closest grid point to desired lat,lon. (query is hardwired for resolution of .3125 deg.) # # Assumption: if we are within .001 of a degree (~80m?), we assume closest gridpoint, otherwise closest 4 are used. # @stmts=("select lon,lat, uncompress(datetime) as datetime, uncompress(prate) as prate,uncompress(tmp2m) as tmp2m, uncompress(q2m) as q2m, uncompress(pressfc) as pressfc, uncompress(uwnd10m) as uwnd10m,uncompress(vwnd10m) as vwnd10m, uncompress(dswsfc) as dswsfc from cfsr${yr} where lon between ${lon}-.3115 and ${lon}+.3115 and lat between (${lat}-.3115) and (${lat}-0.000001);","select lon,lat, uncompress(datetime) as datetime, uncompress(prate) as prate,uncompress(tmp2m) as tmp2m, uncompress(q2m) as q2m, uncompress(pressfc) as pressfc, uncompress(uwnd10m) as uwnd10m,uncompress(vwnd10m) as vwnd10m, uncompress(dswsfc) as dswsfc from cfsr${yr} where lon between ${lon}-.3115 and ${lon}+.312499999 and lat between (${lat}-0.000002) and (${lat}+.3115);"); foreach $stmt (@stmts){ #print "Query:",$stmt,"\n"; $sth = $dbh->prepare($stmt); ## execute it on the database $sth->execute(); my $end = DateTime->now; #print p,$dbh->errstr,"\n"; my $elapsedtime = ($end->subtract_datetime($start))->seconds; #print p,"q1 time(seconds) : $elapsedtime \n"; #print start_ul; # Loop to build an SQL value string by aligning the previously selected variable vectors. # To understand this it might be better to look at the resulting table t${tmptablename} $valuestr_full=""; while($row_ary = $sth->fetchrow_arrayref){ $row_index=0; OUTTER: foreach my $row ( @$row_ary) { @elements=split /,/,$row; if ($#elements < 4) { #print "len:",$#elements,"\n"; next OUTTER; } $el_index=0; foreach $element (@elements){ if($row_index < 1){ @valuestr[$el_index]="@$row_ary[0],@$row_ary[1],\"".$element."\""; ++$el_index; } else { @valuestr[$el_index]="@valuestr[$el_index],$element"; ++$el_index; } } $row_index++; } $valuestr_full=$valuestr_full."(".join('),(',@valuestr)."),"; #print p,$valuestr,"\n"; } chop($valuestr_full); my $end = DateTime->now; my $elapsedtime = ($end->subtract_datetime($start))->seconds; #print p,"l1 time(seconds) : $elapsedtime \n"; $stmt="insert into t${tmptablename} (lon,lat,datetime, prate, tmp2m, q2m, pressfc, uwnd10m, vwnd10m, dswsfc) VALUES ${valuestr_full}"; $sth = $dbh->prepare($stmt); #print p,"Query:",$stmt,"\n"; $sth->execute(); #print p,$dbh->errstr,"\n"; my $end = DateTime->now; my $elapsedtime = ($end->subtract_datetime($start))->seconds; #print p,"q2 time(seconds) : $elapsedtime \n"; } } # # Calculating wind velocity from u,v components and RH from temperature, specific humidity, and pressure # $stmt="update t${tmptablename} set wind=sqrt(pow(uwnd10m,2)+pow(uwnd10m,2)),rh=100*q2m/(.62198*exp(77.3450+0.0057*tmp2m-7235/tmp2m)/(pow(tmp2m,8.2))/(pressfc-exp(77.3450+0.0057*tmp2m-7235/tmp2m)/(pow(tmp2m,8.2))))"; #print p,"Query:",$stmt,"\n"; $sth = $dbh->prepare($stmt); $sth->execute(); # # Query to get surounding grid points, and number of points. # $stmt="select distinct(concat(lat,',',lon)) from t${tmptablename}"; $sth = $dbh->prepare($stmt); $sth->execute(); # # Simple inverse distance weighting power of 2 # http://en.wikipedia.org/wiki/Inverse_distance_weighting # $stncount=0; while($row_ary = $sth->fetchrow_arrayref){ $stncount++; foreach (@$row_ary){ @latlonval[$stncount]=$_; ($glat,$glon)=split /,/,@latlonval[$stncount]; } @gdist[$stncount]=distance($lat,$lon,$glat,$glon,"K"); @wgdist[$stncount]=1/(distance($lat,$lon,$glat,$glon,"K")^$interppow); #@wgdist[$stncount]=1/(distance($lat,$lon,$glat,$glon,"K")); $totgdistance=$totgdistance+@gdist[$stncount]; $totwgdistance=$totwgdistance+@wgdist[$stncount]; #print p,"$glat,$glon distance:".@gdist[$stncount]." ".@wgdist[$stncount]." ".$totwgdistance." ".$stncount; } #print p,"$stncount\n"; # # Creating a new table with the daily summaries necessary for SWAT. # This step could be removed as this is also the result that goes into the final output. # $stmt="create table ta${tmptablename} ENGINE = MEMORY select concat(lat,',',lon) as latlon,date_format(ADDTIME(datetime,\"${tzoffset}:00:00\"),'%Y-%m-%d') as date, max(tmp2m)-273 as tmx,min(tmp2m)-273 as tmn,sum(prate)*3600 as precip, avg(wind) as wind, avg(rh) as avgrh, max(rh)maxrh, min(rh)minrh, sum(dswsfc*3.6/1000.0) as solar from t${tmptablename} group by latlon,date(ADDTIME(datetime,\"${tzoffset}:00:00\")) order by date(ADDTIME(datetime,\"${tzoffset}:00:00\"))"; $sth = $dbh->prepare($stmt); #print p,"Query:",$stmt,"\n"; $sth->execute(); # # Insert weights based on latlon # $stmt="alter table ta${tmptablename} add column weight float not null "; $sth = $dbh->prepare($stmt); $sth->execute(); for ($count=1;$count<=$stncount;$count++){ #print p,$count; #print p,@latlonval[$count]." ".@wgdist[$count]/$totwgdistance; $weight=@wgdist[$count]/$totwgdistance; $stmt="update ta${tmptablename} set weight=$weight where latlon=\"@latlonval[$count]\""; #print p,"Query:",$stmt,"\n"; $sth = $dbh->prepare($stmt); $sth->execute(); } $stmt="create table tb${tmptablename} ENGINE = MEMORY select date, sum(weight*tmx) as tmx,sum(weight*tmn) as tmn,sum(weight*precip) as precip, sum(weight*wind) as wind, sum(weight*avgrh) as avgrh, sum(weight*maxrh) as maxrh, sum(weight*minrh) as minrh, sum(weight*solar) as solar from ta${tmptablename} group by date"; $sth = $dbh->prepare($stmt); #print p,"Query:",$stmt,"\n"; $sth->execute(); # # Select and build table # $stmt="select * from tb${tmptablename} order by date"; #$stmt="select * from t${tmptablename} order by datetime"; $sth = $dbh->prepare($stmt); #print p,"Query:",$stmt,"\n"; $sth->execute(); my $end = DateTime->now; my $elapsedtime = ($end->subtract_datetime($start))->seconds; print p,"Processing time(seconds) : $elapsedtime \n"; $column_names=$sth->{NAME_uc}; print '
$_ | "; } print "|
$_ | "; } else { print "".sprintf("%02.2f",$_ )." | "; } } print "