) {
$r->print($_) unless /^\<\?/;
}
close XML;
}
}
} else {
$r->print($q->header( -type => 'text/html',
-cookie => $cookie ),
$q->start_html('CLIME XML tool'));
my $submit = $q->param('submit');
my $network = $q->param('network');
my $variables = variables($r,$network) if $network;
if ($variables) {
model_editor($r,$q,$ua,$network,$variables);
} else {
mainmenu($r,$q);
}
$r->print($q->end_html);
}
undef $ua;
}
sub xml_header {
my($r,$q) = @_;
$r->print($q->header( -type => 'text/xml' ));
$r->print(''."\n");
}
sub mainmenu {
my($r,$q) = @_;
$r->print(href('root',$me),' ',
href('WMS GetCapabilities',"$me?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetCapabilities"),' ',
href('LakeSet',"$me?doc=LakeSet"),' ',
href('ClimateVarsSet',"$me?doc=ClimateVarsSet"),' ',
href('TopicSet',"$me?doc=TopicSet"),' ',
href('Arrow',"$me?doc=Arrow"),' ',
href('Climate',"$me?doc=Climate"),' ',
href('Patterns',"$me?doc=Patterns"),' ',
"\n");
# dialog box to select a network and proceed
my $networks = networks();
$r->print('',$q->start_form(-method=>'GET'),
$q->popup_menu(-name=>'network',
-values=>$networks),
$q->submit('cmd','Download the model'),
$q->submit('cmd','Document the selected model'),
$q->submit('cmd','Delete the selected model'),
$q->end_form,'
');
}
sub CLIME_DSS {
my($r,$q,$ua) = @_;
xml_element($r,'CLIME_DSS');
xml_element($r,'server',0,'http://clime.tkk.fi/wms.pl?');
xml_element($r,'Location',0,'Europe');
my %help;
my $sth = sql($r,$ua->{dbh},"select name,descr from dsm_helptext");
if ($sth) { for (1..$sth->rows) {
my ($name,$descr) = $sth->fetchrow_array;
$help{$name} = $descr;
}}
for my $k (keys %help) {
xml_element($r,$k,0,$help{$k});
}
xml_element($r,'/CLIME_DSS');
}
sub GetCapabilities_1_1_1 {
my($r,$q) = @_;
xml_element($r,'WMT_MS_Capabilities',
{version=>'1.1.1', updateSequence=>'0'});
xml_element($r,'Service');
xml_element($r,'Name',0,'OGC:WMS');
xml_element($r,'Title',0,'WMS Server for CLIME DSS');
xml_element($r,'Abstract',0,
'Contains maps related to lake management and climate change.');
xml_element($r,'OnlineResource',
{'xmlns:xlink'=>'http://www.w3.org/1999/xlink',
'xlink:type'=>'simple', 'xlink:href'=>$conf->{URLRoot}.$me.'?'},'/');
xml_element($r,'/Service');
xml_element($r,'Capability');
xml_element($r,'Request');
xml_element($r,'GetCapabilities');
xml_element($r,'Format',0,'application/vnd.ogc.wms_xml');
xml_elements($r, ['DCPType'], ['HTTP'], ['Get']);
xml_element($r,'OnlineResource',
{'xmlns:xlink'=>'http://www.w3.org/1999/xlink',
'xlink:type'=>'simple', 'xlink:href'=>$conf->{URLRoot}.$me.'?'},'/');
xml_elements($r, ['/Get'], ['/HTTP'], ['/DCPType'], ['/GetCapabilities']);
xml_elements($r, ['GetMap'], ['Format',0,'image/png'], ['DCPType'], ['HTTP'], ['Get']);
xml_element($r,'OnlineResource',
{'xmlns:xlink'=>'http://www.w3.org/1999/xlink',
'xlink:type'=>'simple', 'xlink:href'=>$conf->{URLRoot}.$me.'?'},'/');
xml_elements($r, ['/Get'], ['/HTTP'], ['/DCPType'], ['/GetMap']);
xml_elements($r, ['GetFeatureInfo'], ['Format',0,'text/html'], ['DCPType'], ['HTTP'], ['Get']);
xml_element($r,'OnlineResource',
{'xmlns:xlink'=>'http://www.w3.org/1999/xlink',
'xlink:type'=>'simple', 'xlink:href'=>$conf->{URLRoot}.$me.'?'},'/');
xml_elements($r, ['/Get'], ['/HTTP'], ['/DCPType'], ['/GetFeatureInfo']);
xml_element($r,'/Request');
xml_element($r,'Exception');
xml_element($r,'Format',0,'application/vnd.ogc.se_blank');
xml_element($r,'/Exception');
xml_element($r,'VendorSpecificCapabilities',0,'/');
xml_elements($r, ['Layer'], ['Name',0,'Europe'], ['Title',0,'Europe']);
# my $g = new Geo::Raster filename=>$conf->{dir}.'/dss/grids/RCAO_g.bil',load=>1;
my ($minX,$minY,$maxX,$maxY) = (800000,700000,7500000,6500000);
my $proj = Geo::Proj4->new(proj => 'laea',
lon_0 => 10,
lat_0 => 52,
ellps => 'GRS80',
x_0 => 4321000,
y_0 => 3210000,
units => 'm');
my($min_lat,$min_lon) = $proj->inverse($minX, $minY);
my($max_lat,$max_lon) = $proj->inverse($maxX, $maxY);
xml_element($r,'SRS',0,'EPSG:3035');
xml_element($r,'LatLonBoundingBox',
{minx=>$min_lon,miny=>$min_lat,maxx=>$max_lon,maxy=>$max_lat},'/');
xml_element($r,'BoundingBox',
{SRS=>'EPSG:3035',
minx=>$minX,miny=>$minY,maxx=>$maxX,maxy=>$maxY},'/');
simple_wms_layer($r,
Name=>'Land cover',
Title=>'Land cover map, CORINE data (c)EEA, Copenhagen, 2005',
Abstract=>'Land cover map based on CORINE data.');
simple_wms_layer($r,
Name=>'Physical',
Title=>'Physical map, includes inland water from CORINE data (c)EEA, Copenhagen, 2005',
Abstract=>'Physical map made mainly from GTOPO30 data.');
xml_elements($r, ['/Layer']);
xml_elements($r, ['Layer'], ['Name',0,'NEAmerica'], ['Title',0,'NEAmerica']);
($minX,$minY,$maxX,$maxY) = (-1644180.228, 4428379.294, 1782695.514, 9997964.943);
$proj = Geo::Proj4->new(proj => 'utm',
zone => 18,
ellps => 'GRS80',
towgs84 => "0,0,0,0,0,0,0",
units => 'm');
($min_lat,$min_lon) = $proj->inverse($minX, $minY);
($max_lat,$max_lon) = $proj->inverse($maxX, $maxY);
xml_element($r,'SRS',0,'EPSG:2149');
xml_element($r,'LatLonBoundingBox',
{minx=>$min_lon,miny=>$min_lat,maxx=>$max_lon,maxy=>$max_lat},'/');
xml_element($r,'BoundingBox',
{SRS=>'EPSG:2149',
minx=>$minX,miny=>$minY,maxx=>$maxX,maxy=>$maxY},'/');
simple_wms_layer($r,
Name=>'NEAmerica Physical',
Title=>'Physical map, mainly from GTOPO30',
Abstract=>'Physical map, mainly from GTOPO30');
xml_elements($r, ['/Layer']);
xml_elements($r, ['/Capability'], ['/WMT_MS_Capabilities']);
}
sub simple_wms_layer {
my $r = shift;
my %p = @_;
xml_element($r,'Layer');
xml_element($r,'Name',0,$p{Name});
xml_element($r,'Title',0,$p{Title});
xml_element($r,'Abstract',0,$p{Abstract});
xml_element($r,'/Layer');
}
sub GetMap_1_1_1 {
my($r,$q) = @_;
my $layers = my_param($q, 'LAYERS');
$layers = my_param($q, 'Layers') unless $layers;
my @layers = split/,/, $layers;
my $srs = my_param($q,'SRS'); # we serve only 3035
my $epsg;
if ($srs =~ /EPSG:(\d+)/) {
$epsg = $1;
}
my @bbox1 = split/,/, my_param($q,'BBOX');
my @bbox;
for (@bbox1) {
my($a) = $_ =~ /([\d\.\-E]+)/;
push @bbox,$a;
}
my ($w) = my_param($q,'WIDTH') =~ /(\d+)/;
my ($h) = my_param($q,'HEIGHT') =~ /(\d+)/;
my $format = my_param($q,'FORMAT');
$w = 523 unless $w;
$h = 453 unless $h;
# jotain rajaa:
$w = 2000 if $w > 2000;
$h = 2000 if $h > 2000;
my $name = time();
my $clip = '/l/bin/gdal_translate -of PNG -not_strict';
my $raster = "$conf->{dir}/dss/CORINE/lceugr100_00_pct.tif";
my $vector = "$conf->{dir}/dss/CORINE/rajat";
my $shoreline = "$conf->{dir}/dss/CORINE/gshhs_land";
my $png = "$conf->{dir}/dss/tmp/$name.png";
my @bgcolor = (138,200,199); # sea
# system "$clip -outsize $w $h -projwin $bbox[0] $bbox[3] $bbox[2] $bbox[1] $raster $png";
my $pixel_size = ($bbox[2] - $bbox[0])/$w; # assuming..
my ($minX,$minY,$maxX,$maxY) = @bbox;
@bbox = ($minX,$minY,$minX+$w*$pixel_size,$minY+$h*$pixel_size);
my @wms = ();
#print STDERR "$layers\n";
if ($layers =~ /^Land/) {
# country fills
$wms[0] = new Geo::Vector $vector unless defined $wms[0];
$wms[0]->{alpha} = 255;
$wms[0]->{COLOR} = [255,255,255,255]; # no data
# CORINE
$wms[1] = new Geo::Raster $raster unless defined $wms[1];
$wms[1]->{GDAL}->{dataset}->GetRasterBand(1)->SetNoDataValue(255);
# country borders
$wms[2] = new Geo::Vector $vector;
# $wms[2] = new Geo::Vector $shoreline;
$wms[2]->{RenderAs} = 2;
$wms[2]->{alpha} = 255;
$wms[2]->{COLOR} = [0,0,0,255];
# catchments
$wms[3] = new Geo::Vector "$conf->{dir}/dss/CORINE/lakes_t";
$wms[3]->{alpha} = 100;
$wms[3]->{COLOR} = [100,255,100,255];
$wms[4] = new Geo::Vector "$conf->{dir}/dss/CORINE/lakes_t";
$wms[4]->{alpha} = 130;
$wms[4]->{COLOR} = [0,0,0,255];
$wms[4]->{RenderAs} = 2;
# print STDERR "$bbox[0], $bbox[2] ",$bbox[2] - $bbox[0]," $w\n";
# print STDERR "$bbox[0], $bbox[3] $pixel_size\n";
} elsif ($layers =~ /^Physical/) {
$wms[0] = new Geo::Raster "$conf->{dir}/dss/CORINE/elev.tiff";
my $corine = new Geo::Raster $raster;
$corine->cache(@bbox,$pixel_size);
$wms[1] = $corine == 41; # lakes
$wms[1]->nodata_value(0);
my $ct = $wms[1]->get_color_table(1);
my @c = (21,40,237,255);
#my @c = (128,242,230,255); # same as corine 41
$ct->SetColorEntry(1,\@c);
$wms[2] = new Geo::Vector $vector;
$wms[2]->{RenderAs} = 2;
$wms[2]->{alpha} = 255;
$wms[2]->{COLOR} = [0,0,0,255];
$wms[3] = new Geo::Vector "$conf->{dir}/dss/CORINE/lakes_t";
$wms[3]->{alpha} = 100;
$wms[3]->{COLOR} = [100,255,100,255];
$wms[4] = new Geo::Vector "$conf->{dir}/dss/CORINE/lakes_t";
$wms[4]->{alpha} = 130;
$wms[4]->{COLOR} = [0,0,0,255];
$wms[4]->{RenderAs} = 2;
} else { # NE America
$wms[0] = new Geo::Raster "$conf->{dir}/dss/NEAmerica/elev.tiff";
$wms[1] = new Geo::Vector "$conf->{dir}/dss/NEAmerica/PAL";
$wms[1]->{RenderAs} = 2;
$wms[1]->{alpha} = 200;
$wms[1]->{COLOR} = [255,0,0,255];
$wms[2] = new Geo::Vector "$conf->{dir}/dss/NEAmerica/NHDWaterbody";
$wms[2]->{alpha} = 200;
$wms[2]->{COLOR} = [21,40,237,255];
}
my $pixbuf = Gtk2::Ex::Geo::Renderer->new
([@wms],$bbox[0],$bbox[3],$pixel_size,$w,$h,
0,0,@bgcolor);
$pixbuf->save($png,'png');
serve_document($r,$q,$mime_types,$png);
}
sub my_param {
my $p = $_[0]->param($_[1]);
return defined $p ? $p : '';
}
sub xml_elements {
my $r = shift;
for (@_) {
xml_element($r,@$_);
}
}
sub xml_element {
my($r,$element,$attributes,$content) = @_;
$r->print("<$element");
if ($attributes) {
for my $a (keys %$attributes) {
$r->print(" $a=\"$attributes->{$a}\"");
}
}
unless (defined $content) {
$r->print(">\n");
} elsif ($content eq '/') {
$r->print("/>\n");
} else {
$r->print(">$content$element>\n");
}
}
sub model_editor {
my($r,$q,$ua,$network,$variables) = @_;
my $distributions = distributions();
my $mapping = mapping($r,$ua->{dbh},$network) if $network;
my @topics;
my @csites;
my $sth = sql($r,$ua->{dbh},"select name from dsm_water_quality_issue order by name");
if ($sth) { for (1..$sth->rows) {
my ($name) = $sth->fetchrow_array;
push @topics,$name;
}}
$sth = sql($r,$ua->{dbh},"select lo.name from dsm_lake l join dsm_location lo on l.id=lo.id order by lo.name");
if ($sth) { for (1..$sth->rows) {
my ($name) = $sth->fetchrow_array;
push @csites,$name;
}}
if ($network and $q->param('cmd') =~ /^Store this/) {
my $new_mapping = { topic => $q->param("topic"),
csite => $q->param("csite"),
descr => $q->param("descr"),
};
for my $i (0..$#$variables) {
my $var = $variables->[$i];
$new_mapping->{$var}{sp_dist} = $q->param("var $i");
$new_mapping->{$var}{descr} = $q->param("descr var $i");
}
store_mapping($r,$q,$ua->{dbh},$mapping,$new_mapping,$network);
$r->print("The network documentation was updated.
\n");
mainmenu($r,$q);
} elsif ($network and $q->param('cmd') =~ /^Delete the selected model/) {
my ($filename) = $network =~ /^([\w\.]+)/;
rename "$conf->{dir}/dss/networks/$filename","$conf->{dir}/dss/networks/$filename.saved";
sql($r,$ua->{dbh},"delete from bn_map where network='$network'");
mainmenu($r,$q);
} elsif ($q->param('cmd') =~ /^Document the selected model/) {
# now ask map variable to distribution
unshift @$distributions,'Nothing';
$r->print("Now documenting network '$network' (",href('back',$me),")
\n",'',
$q->start_form,
$q->hidden('network',$network));
my $descr = exists $mapping->{descr} ? $mapping->{descr} : '';
$r->print("Description of $network:
\n",
$q->textarea(-name=>"descr",
-default=>$descr,
-rows=>5,
-columns=>40),"
\n");
my $topic = exists $mapping->{topic} ? $mapping->{topic} : 'Nothing';
$r->print("The topic of this network is ",
$q->popup_menu(-name=>"topic",
-values=>\@topics,
-default=>$topic),"
\n");
my $csite = exists $mapping->{csite} ? $mapping->{csite} : 'Nothing';
$r->print("The calibration site of this network is ",
$q->popup_menu(-name=>"csite",
-values=>\@csites,
-default=>$csite),"
\n");
for my $i (0..$#$variables) {
my $variable = $variables->[$i];
my $default = exists $mapping->{$variable} ? $mapping->{$variable}{sp_dist} : 'Nothing';
my $descr = exists $mapping->{$variable} ? $mapping->{$variable}{descr} : '';
$r->print("Description of $variable:
\n",
$q->textarea(-name=>"descr var $i",
-default=>$descr,
-rows=>5,
-columns=>40),"
\n");
$r->print("Link $variable to ",
$q->popup_menu(-name=>"var $i",
-values=>$distributions,
-default=>$default),"
\n");
}
$r->print('
',$q->submit('cmd','Store this'),'
',$q->end_form);
}
}
sub networks {
my @networks;
if (opendir(DIR, "$conf->{dir}/dss/networks/")) {
@networks = grep { /\.xml$/i && -r "$conf->{dir}/dss/networks/$_" } sort readdir(DIR);
closedir DIR;
}
return \@networks;
}
sub variables {
my($r,$network,$only_inputs) = @_;
# SHOULD start using XML::LibXML
my $parser = new XML::DOM::Parser;
my $doc;
eval {
$doc = $parser->parsefile($conf->{dir}.'/dss/networks/'.$network);
};
if ($@ or !$doc) {
$r->print("There was an error in parsing $network:
$@
");
return 0;
} else {
my %variables;
my @variables = $doc->getElementsByTagName('VARIABLE');
for my $variable (@variables) {
my @list = $variable->getChildNodes;
for my $n (@list) {
my $name = $n->getNodeName;
if ($name eq 'NAME') {
my $v = $n->getFirstChild()->getNodeValue;
$variables{$v} = 1;
}
}
}
if ($only_inputs) {
@variables = $doc->getElementsByTagName('DEFINITION');
for my $variable (@variables) {
my @list = $variable->getChildNodes;
my $for;
my $given;
for my $n (@list) {
my $name = $n->getNodeName;
if ($name eq 'FOR') {
$for = $n->getFirstChild()->getNodeValue;
} elsif ($name eq 'GIVEN') {
$given = $n->getFirstChild()->getNodeValue;
}
}
if ($for and $given) {
delete $variables{$for};
}
}
}
@variables = sort keys %variables;
return \@variables;
}
}
sub distributions {
my($full_name) = @_;
my @distributions;
if (opendir(DIR, "$conf->{dir}/dss/distributions")) {
my @x = grep { /\.png$/i && -r "$conf->{dir}/dss/distributions/$_" } readdir(DIR);
my %x;
for (@x) {
unless ($full_name) {
$_ =~ s/(RCAO-H|RCAO-E|HadRM3p)_(C|A2|B2)_//;
$_ =~ s/\.png$//i;
}
$x{$_} = 1;
}
@distributions = sort keys %x;
closedir DIR;
}
return \@distributions;
}
sub store_mapping {
my($r,$q,$dbh,$old_mapping,$new_mapping,$network) = @_;
for my $var (keys %$new_mapping) {
my ($sp_dist, $descr) = ('','');
if ($var eq 'topic' or $var eq 'csite') {
$sp_dist = $new_mapping->{$var};
} elsif ($var eq 'descr') {
$descr = $new_mapping->{descr};
} else {
$sp_dist = $new_mapping->{$var}{sp_dist};
$descr = $new_mapping->{$var}{descr};
}
if (exists $old_mapping->{$var}) {
sql($r,$ua->{dbh},"update bn_map set sp_dist='$sp_dist',descr='$descr' ".
"where network='$network' and varname='$var'");
} else {
sql($r,$ua->{dbh},"insert into bn_map (network,varname,sp_dist,descr) ".
"values ('$network','$var','$sp_dist','$descr')");
}
}
}
sub mapping {
my($r,$dbh,$network) = @_;
my $sth = sql($r,$dbh,"select varname,sp_dist,descr from bn_map where network='$network'");
my $mapping = {};
if ($sth) { for (1..$sth->rows) {
my ($var,$sp_dist,$descr) = $sth->fetchrow_array;
if ($var eq 'topic' or $var eq 'csite') {
$mapping->{$var} = $sp_dist;
} elsif ($var eq 'descr') {
$mapping->{$var} = $descr;
} else {
$mapping->{$var}{sp_dist} = $sp_dist;
$mapping->{$var}{descr} = $descr;
}
}}
return $mapping;
}
sub mappings {
my($r,$dbh) = @_;
my $sth = sql($r,$dbh,"select network,varname,sp_dist,descr from bn_map");
my $mappings = {};
if ($sth) { for (1..$sth->rows) {
my ($network,$var,$sp_dist,$descr) = $sth->fetchrow_array;
if ($var eq 'topic' or $var eq 'csite') {
$mappings->{$network}{$var} = $sp_dist;
} elsif ($var eq 'descr') {
$mappings->{$network}{$var} = $descr;
} else {
$mappings->{$network}{$var}{sp_dist} = $sp_dist;
$mappings->{$network}{$var}{descr} = $descr;
}
}}
return $mappings;
}
sub TopicSet {
my($r,$q,$dbh,$u) = @_;
$r->print("");
my $mappings = mappings($r,$dbh);
my $distributions = 0; #distributions(1);
my $TopicCode = 0;
my %grid;
$grid{'RCAO'} = new Geo::Raster filename=>$conf->{dir}.'/dss/grids/RCAO_g.bil',load=>1;
$grid{'HadRM3p'} = new Geo::Raster filename=>$conf->{dir}.'/dss/grids/HadRM3p_s_g.bil',load=>1;
my $proj = Geo::Proj4->new(proj => 'laea',
lon_0 => 10,
lat_0 => 52,
ellps => 'GRS80',
x_0 => 4321000,
y_0 => 3210000,
units => 'm');
my %rlatlon;
my $sql =
'select dsm_location.id,dsm_location.name,dsm_lake.longitude,dsm_lake.latitude ' .
'from dsm_lake, dsm_location where dsm_lake.id=dsm_location.id and dsm_lake.latitude > 0 order by dsm_location.name';
my $sth = sql($r,$dbh,$sql);
if ($sth) {
for (1..$sth->rows) {
my ($id,$lake,$lon,$lat) = $sth->fetchrow_array;
next unless $lon and $lat;
my @xy = $proj->forward($lat,$lon);
for my $model (keys %grid) {
my ($minX,$minY,$maxX,$maxY) = $grid{$model}->world();
next if $xy[0] < $minX;
next if $xy[1] < $minY;
next if $xy[0] > $maxX;
next if $xy[1] > $maxY;
my @ij = $grid{$model}->w2g(@xy);
my $r = $grid{$model}->get(@ij);
next unless defined $r;
my $rlat = int($r/100);
my $rlon = $r % 100;
$rlat = 0 unless $rlat;
$rlon = 0 unless $rlon;
# print STDERR "$lake,$model,$r,$rlat,$rlon,@ij\n";
$rlatlon{$lake}{$model}{rlat} = $rlat;
$rlatlon{$lake}{$model}{rlon} = $rlon;
}
}
}
my %topics;
$sql = 'select name,descr from dsm_water_quality_issue';
$sth = sql($r,$dbh,$sql);
if ($sth) {
for (1..$sth->rows) {
my ($name,$descr) = $sth->fetchrow_array;
$topics{$name} = $descr;
}
}
for my $network (sort keys %$mappings) {
$TopicCode++;
$r->print("");
my $n = $network;
$n =~ s/\.xml$//i;
my @n = split /_/,$n;
if ($mappings->{$network}{topic}) {
$r->print("$mappings->{$network}{topic}");
} else {
$r->print("$n[0]");
}
$r->print("$n[1]");
$r->print("$mappings->{$network}{descr}") if $mappings->{$network}{descr};
$r->print("$network");
for my $model ('RCAO-E','RCAO-H','HadRM3p') {
my $g = $model;
$g =~ s/\-.*$//;
my $lake = $mappings->{$network}{csite};
$r->print("") unless $lake;
next unless $lake;
my $m = substr($model,0,1).substr($model,-1,1);
my $fn = $m."_".$rlatlon{$lake}{$g}{rlat}."_".$rlatlon{$lake}{$g}{rlon}.".png";
my $file = $conf->{dir}.'/dss/sims/'.$fn;
if (-r $file) {
$r->print("",
"$conf->{URLRoot}$me?SERVICE=file\&\;REQUEST=sims\&\;doc=$fn",
"");
} else {
$r->print("$file");
}
}
# $r->print("$conf->{URLRoot}?dss=sims/RE_6065.png");
# $r->print("$conf->{URLRoot}?dss=sims/RH_6065.png");
# $r->print("$conf->{URLRoot}?dss=networks/$network");
$r->print("$conf->{URLRoot}$me?network=$network\&\;cmd=Download+the+model");
my $parser = new XML::DOM::Parser;
my $doc;
eval {
$doc = $parser->parsefile($conf->{dir}.'/dss/networks/'.$network);
};
if ($@ or !$doc) {
$r->print("$network$@");
$r->print("");
next;
}
my @ClimateScenarios;
my $sql =
'select name ' .
'from dsm_cc_scenario order by name';
my $sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($scenarioName) = $sth->fetchrow_array;
push(@ClimateScenarios,$scenarioName);
}
}
my @ClimateModels;
$sql =
'select name ' .
'from dsm_model where system = 4234 order by name';
$sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($ModelName) = $sth->fetchrow_array;
push(@ClimateModels,$ModelName);
}
}
my @variables = $doc->getElementsByTagName('VARIABLE');
for my $variable (@variables) {
my @list = $variable->getChildNodes;
for my $n (@list) {
my $name = $n->getNodeName;
if ($name eq 'NAME') {
my $v = $n->getFirstChild()->getNodeValue;
$r->print("");
my $vdata = $mappings->{$network}{$v};
my $descr = $vdata->{descr} if $vdata;
$r->print("$descr");
my $sp_dist = $vdata->{sp_dist} if $vdata;
unless ($sp_dist) {
$r->print("");
next;
}
$r->print("$sp_dist");
my ($prefix);
if ($sp_dist =~ /^\w/) {
($prefix) = $sp_dist =~ /^([\d\-a-zA-Z]+_)/;
$prefix = '' unless $prefix;
$sp_dist =~ s/^([\d\-a-zA-Z]+_)//;
} else {
$prefix = '';
}
$r->print("$sp_dist");
$r->print("$prefix");
my $count = 0;
if (0) { # not used any more
for my $fn (@$distributions) {
# $r->print("$fn");
for my $model (@ClimateModels) {
for my $scenario (@ClimateScenarios) {
# if ($fn =~ /$sp_dist/) {
if ($fn eq $prefix.$model.'_'.$scenario.'_'.$sp_dist.'.png') {
$r->print("",
"$conf->{URLRoot}?dss=distributions/$fn");
$count++;
}
}
}
}
$r->print("") unless $count;
}
$r->print(" ");
}
}
}
$doc->dispose;
$r->print("");
}
$r->print("");
}
sub ClimateVarsSet {
my($r,$q,$dbh,$u) = @_;
$r->print("");
my ($minX,$minY,$maxX,$maxY,$EPSG);
my $loc = $q->param('Location');
if ($loc eq 'NEAmerica') {
($minX,$minY,$maxX,$maxY) = (-1644180.228, 4428379.294, 1782695.514, 9997964.943);
$EPSG = 2149;
} else {
my $grid = new Geo::Raster filename=>$conf->{dir}.'/dss/grids/RCAO_g.bil',load=>1;
($minX,$minY,$maxX,$maxY) = $grid->world();
$EPSG = 3035;
}
xml_element($r,'BoundingBox',
{SRS=>'EPSG:'.$EPSG,
minx=>$minX,miny=>$minY,maxx=>$maxX,maxy=>$maxY},'/');
my %ext = (1=>'png',2=>'jpg');
my %subdir = (1 => 'distributions', 2 => 'dbv', 3=> 'rdbv');
my @ClimateScenarios;
my %scenarioHelp;
my $sql =
'select name, descr ' .
'from dsm_cc_scenario order by name';
my $sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($name,$descr) = $sth->fetchrow_array;
push(@ClimateScenarios,$name);
$scenarioHelp{$name} = $descr;
}
}
my @ClimateModels;
my %modelsHelp;
$sql =
'select name, descr ' .
'from dsm_model where system = 4234 order by name';
$sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($name,$descr) = $sth->fetchrow_array;
push(@ClimateModels,$name);
$modelsHelp{$name} = $descr;
}
}
my @variables;
my %min;
my %max;
my %scale;
my %unit;
my %variableHelp;
my %fname;
my %rules;
$sql =
'select v.acronym, v.name, v.descr, v.scale, u.name, vc.rule, c.acronym '.
'from dsm_unit u join dsm_variable v on u.id=v.unit '.
'left join dsm_variable_category vc on v.id=vc.var '.
'left join dsm_category c on c.id=vc.cat';
$sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($acro,$name,$descr,$scale,$unit,$rule,$cat) = $sth->fetchrow_array;
$scale{$acro} = $scale;
$fname{$acro} = $name;
$unit{$acro} = $unit;
$variableHelp{$acro} = $descr;
push @{$rules{$acro}},{$rule=>$cat};
}
}
@variables = sort keys %scale;
my @seasons;
my %seasonHelp;
$sql = 'select name,descr from dsm_season order by number';
$sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($name,$descr) = $sth->fetchrow_array;
push(@seasons,$name);
$seasonHelp{$name} = $descr;
}
}
my %ClimateScenarioID;
$r->print("");
my $id = 1;
for my $scenario (@ClimateScenarios) {
# next if $scenario eq 'C'; # needed!
$r->print("");
$ClimateScenarioID{$scenario} = $id;
$id++;
$r->print("$scenario");
$r->print("$scenarioHelp{$scenario}");
$r->print("");
}
$r->print("");
my %ClimateModelID;
$r->print("");
$id = 1;
for my $model (@ClimateModels) {
$r->print("");
$ClimateModelID{$model} = $id;
$id++;
$r->print("$model");
my $grid = $model =~ /^Had/ ? 'HadRM3p_s' : 'RCAO';
$r->print("$conf->{URLRoot}$me?SERVICE=file\&\;REQUEST=grids\&\;doc=${grid}_g.png");
$r->print("$modelsHelp{$model}");
$r->print("");
}
$r->print("\n");
my %SeasonID;
$r->print("");
$id = 1;
for my $season (@seasons) {
$r->print("");
$SeasonID{$season} = $id;
$id++;
$r->print("$season");
$r->print("$seasonHelp{$season}");
$r->print("");
}
$r->print("\n");
my @MapTypes = ("Absolute change", "Relative change");
my %mapTypeHelp;
$mapTypeHelp{"Absolute change"} = "The absolute change map represents the change of the selected climate variable in absolute values.";
$mapTypeHelp{"Relative change"} = "The relative change map represents the change of the selected climate variable as a percentage.";
my %MapTypeID;
$r->print("");
$id = 1;
for my $mapType (@MapTypes) {
$r->print("");
$MapTypeID{$mapType} = $id;
$id++;
$r->print("$mapType");
$r->print("$mapTypeHelp{$mapType}");
$r->print("");
}
$r->print("\n");
my %VariableID;
$r->print("");
$id = 1;
for my $variable (@variables) {
$r->print("");
$VariableID{$variable} = $id;
$id++;
$r->print("$variable\n");
$r->print("$scale{$variable}\n");
$r->print("$unit{$variable}\n");
$r->print("$fname{$variable}\n");
$r->print("$variableHelp{$variable}\n");
for my $rule (@{$rules{$variable}}) {
for my $key (keys %$rule) {
my $k = $key;
$k =~ s/\</g;
$k =~ s/\>/>/g;
$r->print("{$key}\">$k\n");
}
}
$r->print("\n");
}
$r->print("\n");
my $URLRoot = $conf->{URLRoot};
for my $variable (@variables) {
$r->print('');
for my $scenario (@ClimateScenarios) {
for my $model (@ClimateModels) {
for my $season (@seasons) {
my $fn = $model.'_'.$scenario.'_'.lc($season).'_'.$variable.'.'.$ext{1};
$r->print("",
$URLRoot,$me,"?SERVICE=file\&\;REQUEST=",$subdir{2},"\&\;doc=",$fn,
"");
}
for my $season (@seasons) {
my $fn = $model.'_'.$scenario.'_'.lc($season).'_'.$variable.'.'.$ext{1};
$r->print("",
$URLRoot,$me,"?SERVICE=file\&\;REQUEST=",$subdir{3},"\&\;doc=",$fn,
"");
}
}
}
$r->print("");
}
$r->print("");
}
sub MapSet {
my($r,$q,$dbh,$u) = @_;
$r->print(''."\n");
$r->print(''."\n");
$r->print(''."\n");
$r->print(''."\n");
$r->print(''."\n");
}
sub LakeData {
my($r,$q,$dbh,$u) = @_;
my ($lake_id) = my_param($q,'dbid') =~ /(\d+)/;
unless ($lake_id) {
xml_element($r,'Error',0,'lake db id missing');
return;
}
my $name;
my $sql = 'select l.id,l.name,s.longitude,s.latitude ' .
"from dsm_lake s, dsm_location l where s.id=l.id and s.latitude > 0 and s.id='$lake_id'";
my $sth = sql($r,$dbh,$sql);
if ($sth) { for (1..$sth->rows) {
my ($id,$lake,$lon,$lat) = $sth->fetchrow_array;
$name = $lake;
}}
unless ($name) {
xml_element($r,'Error',0,'no such lake');
return;
}
my $catchment;
my $boundary;
my $catchments = new Geo::Vector("$conf->{dir}/dss/CORINE/lakes_t");
$catchments->{ogr_layer}->ResetReading();
my $test = $name;
$test =~ s/[^a-zA-Z]//g;
my $feature;
while ($feature = $catchments->{ogr_layer}->GetNextFeature) {
my $n = $feature->GetFieldAsString('Name');
$n =~ s/[^a-zA-Z]//g;
if ($n eq $test) {
$catchment = $feature->GetFID();
my $geom = $feature->GetGeometryRef();
$boundary = $geom->GetBoundary;
}
}
unless (defined $catchment) {
xml_element($r,'Error',0,'no catchment for the lake');
return;
}
my $lc = new Geo::Raster "$conf->{dir}/dss/CORINE/lceugr100_00_pct.tif";
$lc->cache($catchments->world(feature=>$catchment));
unless ($lc->{M} and $lc->{N}) {
xml_element($r,'Error',0,'no non-null catchment for the lake');
return;
}
my $clip = $catchments->rasterize(like=>$lc, feature=>$catchment);
$lc *= $clip;
# mapping to land cover classes
my %map;
my %classes = (301=>'Artificial surface',
302=>'Agricultural areas',
303=>'Peat contributing areas',
304=>'Forests',
305=>'Wetlands',
306=>'Inland waters',
307=>'Marine waters');
for my $i (1..11) {$map{$i} = 301}
for my $i (12..22) {$map{$i} = 302}
for my $i (24..27,29,35,36) {$map{$i} = 303}
for my $i (23,28,30..34) {$map{$i} = 304}
for my $i (37..39) {$map{$i} = 305}
for my $i (40..41) {$map{$i} = 306}
for my $i (42..44) {$map{$i} = 307}
$lc->map(\%map);
my $contents = $lc->contents;
my $cell_area = 1; # ha
xml_element($r,'LandCover');
xml_element($r,'ForCatchmentOfLake',0,$name);
my $sum = 0;
for my $v (keys %$contents) {
next if $v == 0;
my $a = $cell_area*$contents->{$v};
$sum += $cell_area*$contents->{$v};
}
for my $v (sort {$contents->{$b}<=>$contents->{$a}} keys %$contents) {
next if $v == 0;
my $unit = 'ha';
my $p = $cell_area*$contents->{$v}/$sum*100;
# my $a = $cell_area*$contents->{$v};
# if ($a > 100) {
# $a /= 100;
# $unit = 'km2';
# }
my $c = $classes{$v};
$c = 'not classified' unless $c;
#xml_element($r,'Amount',{LandCoverClass=>$c,Unit=>$unit},$a);
xml_element($r,'Amount',{LandCoverClass=>$c,Unit=>'%'},$p);
}
my $unit = 'ha';
if ($sum > 100) {
$sum /= 100;
$unit = 'km2';
}
xml_element($r,'Total',{Unit=>$unit},$sum);
xml_element($r,'/LandCover');
}
sub LakeSet {
my($r,$q,$dbh,$u) = @_;
my $raster = new Geo::Raster filename=>"$conf->{dir}/dss/lc_clip.bil",load=>1;
my $catchments = my_param($q,'catchments');
my %catchments;
my $layer;
if ($catchments eq 'yes') {
$catchments = ogr::Open("$conf->{dir}/dss/CORINE/");
$layer = $catchments->GetLayerByName('lakes_t');
$layer->ResetReading();
my $feature;
while ($feature = $layer->GetNextFeature) {
my $name = $feature->GetFieldAsString('Name');
$name =~ s/[^a-zA-Z]//g;
$catchments{$name} = $feature->GetFID();
}
}
# EPSG 3035:
my $proj = Geo::Proj4->new(proj => 'laea',
lon_0 => 10,
lat_0 => 52,
ellps => 'GRS80',
x_0 => 4321000,
y_0 => 3210000,
units => 'm');
my $sql = 'select l.id,l.name,s.longitude,s.latitude ' .
'from dsm_lake s, dsm_location l where s.id=l.id and s.latitude > 0 order by l.name';
my $sth = sql($r,$dbh,$sql);
$r->print(''."\n");
my ($minX,$minY,$maxX,$maxY) = $raster->world();
my $sitecode = 1;
if ($sth) { for (1..$sth->rows) {
my ($id,$lake,$lon,$lat) = $sth->fetchrow_array;
next unless $lon and $lat;
my @xy = $proj->forward($lat,$lon);
next if $xy[0] < $minX;
next if $xy[1] < $minY;
next if $xy[0] > $maxX;
next if $xy[1] > $maxY;
# my @ij = $raster->w2g(@xy);
$r->print("\n");
$r->print("$id\n");
$r->print("$lake\n");
# $r->print("$ij[1]$ij[0]\n");
xml_element($r, 'Location');
xml_element($r, 'X', 0, $xy[0]);
xml_element($r, 'Y', 0, $xy[1]);
xml_element($r, '/Location');
xml_element($r, 'Catchment');
my $name = $lake;
$name =~ s/[^a-zA-Z]//g;
if ($layer and $catchments{$name}) {
my $f = $layer->GetFeature($catchments{$name});
my $g = $f->GetGeometryRef();
$g = $g->GetGeometryRef(0);
for my $i (0..$g->GetPointCount-1) {
xml_element($r, 'Point');
xml_element($r, 'X', 0, $g->GetX($i));
xml_element($r, 'Y', 0, $g->GetY($i));
xml_element($r, '/Point');
}
}
xml_element($r, '/Catchment');
$r->print(''."\n");
$sitecode++;
}}
$r->print(''."\n");
}
sub Arrow {
my($r,$q,$dbh,$u) = @_;
my $model = $q->param('model');
my $scenario = $q->param('scenario');
my $rlat = $q->param('rlat');
my $rlon = $q->param('rlon');
$model =~ s/[^\w-]//g;
$scenario =~ s/\W//g;
$rlat =~ s/\W//g;
$rlon =~ s/\W//g;
$r->print('parameters missing'."\n"), return
unless $model and $scenario and $rlat and $rlon;
my $gm = $model;
$gm =~ s/-\w$//;
# $gm .= '_s' if $gm =~ /^H/;
my $sql =
"select a.from_rlat,a.from_rlon,a.sim,g_from.i,g_from.j,g_to.i,g_to.j ".
"from arrows a, grid g_from, grid g_to ".
"where a.model='$model' and a.sc='$scenario' and ".
"a.to_rlat='$rlat' and a.to_rlon='$rlon' and ".
"g_from.model='$gm' and g_from.rlat=a.from_rlat and g_from.rlon=a.from_rlon and ".
"g_to.model='$gm' and g_to.rlat='$rlat' and g_to.rlon='$rlon'";
my $sth = sql($r,$dbh,$sql);
my $from_rlat = '';
my $from_rlon = '';
my ($from_i,$from_j,$to_i,$to_j) = ('','','','');
my $sim = 0;
if ($sth) {
for (1..$sth->rows) {
($from_rlat,$from_rlon,$sim,$from_i,$from_j,$to_i,$to_j) = $sth->fetchrow_array;
}
}
$r->print('');
$r->print("$sql");
$r->print("$model");
$r->print("$scenario");
$r->print("$rlat");
$r->print("$rlon");
$r->print("$from_rlat");
$r->print("$from_rlon");
$r->print("$from_i");
$r->print("$from_j");
$r->print("$to_i");
$r->print("$to_j");
$r->print("$sim");
$r->print(''."\n");
}
sub Climate {
my($r,$q,$dbh,$u) = @_;
my $model = $q->param('model');
my $season = $q->param('season');
my $rlat = $q->param('rlat');
my $rlon = $q->param('rlon');
$model =~ s/[^\w-]//g;
$season =~ s/\W//g;
$rlat =~ s/\W//g;
$rlon =~ s/\W//g;
$r->print('parameters missing'."\n"), return
unless $model and $season and $rlat and $rlon;
my $sql = "select sc,var,v from climate ".
"where model='$model' and season='".lc($season)."' and rlat='$rlat' and rlon='$rlon' ".
"order by sc,var";
# $r->print($sql);
my $sth = sql($r,$dbh,$sql);
$r->print('');
$r->print(" $model ");
$r->print(" $season ");
$r->print(" $rlat ");
$r->print(" $rlon ");
# $r->print("");
my $sc0 = '';
if ($sth) {
for (1..$sth->rows) {
my ($sc,$var,$v) = $sth->fetchrow_array;
if ($sc ne $sc0) {
if ($sc0) {
$r->print("no data") if $model eq 'HadRM3p';
$r->print('');
}
$r->print("");
$sc0 = $sc;
}
$r->print(" $v");
}
$r->print("no data") if $model eq 'HadRM3p';
$r->print('') if $sth->rows;
}
$r->print('');
}
sub Patterns {
my($r,$q,$dbh,$u) = @_;
my $model = $q->param('model');
my $season = $q->param('season');
my $rlat = $q->param('rlat');
my $rlon = $q->param('rlon');
$model =~ s/[^\w-]//g;
$season =~ s/\W//g;
$rlat =~ s/\W//g;
$rlon =~ s/\W//g;
$r->print('parameters missing'."\n"), return
unless $model and $season and $rlat and $rlon;
my $sql = "select sc,cat,n from patterns ".
"where model='$model' and season='".lc($season)."' and rlat='$rlat' and rlon='$rlon' ";
my %data;
my %cats;
my $sth = sql($r,$dbh,$sql);
my @all_variables = ('precip','runoff','snow','soilw','t2m','w10m');
my @variables = ('precip','snow','t2m','w10m');
if ($sth) {
for (1..$sth->rows) {
my ($sc,$cat,$n) = $sth->fetchrow_array;
my %pat;
my @cats = split /,/,$cat;
for my $j (0..$#all_variables) {
my $c = $cats[$j];
$pat{$all_variables[$j]} = $c;
}
@cats = ();
for my $j (0..$#variables) {
push @cats,$pat{$variables[$j]};
}
$cat = join(',',@cats);
$data{$cat}{$sc} += $n;
$cats{$cat} = 0 unless defined $data{$cat};
$cats{$cat} += $n;
}
}
my @patterns = sort {$cats{$b} <=> $cats{$a}} keys %cats;
@patterns = @patterns[0..3];
$r->print('');
#$r->print(" $sql ");
$r->print(" $model ");
$r->print(" $season ");
$r->print(" $rlat ");
$r->print(" $rlon ");
my @scenarios = ('C','A2','B2');
my $i = 1;
for my $pattern (@patterns) {
$r->print("");
my @cats = split /,/,$pattern;
for my $j (0..$#variables) {
my $c = $cats[$j];
$c = 'yes' if $c eq 'snow';
$c = 'no' if $c eq 'no_snow';
$c = 'no data' if $variables[$j] eq 'snow' and $model eq 'HadRM3p';
$r->print("<$variables[$j]>$c$variables[$j]>");
}
for my $scenario (@scenarios) {
my $p = $data{$pattern}{$scenario};
$p = 0 unless defined $p;
$p = sprintf("%.0f",$p/90*100) . '%';
$r->print("<$scenario>$p$scenario>");
}
$i++;
$r->print("");
}
$r->print(''."\n");
}
sub ClimateTimeseries {
my($r,$q,$dbh,$u) = @_;
my $model = $q->param('model');
my $variable = $q->param('variable');
my $scenario = $q->param('scenario');
my $rlat = $q->param('rlat');
my $rlon = $q->param('rlon');
$model =~ s/[^\w-]//g;
$variable =~ s/[^\w\s-]//g;
$scenario =~ s/[^\w-]//g;
$rlat =~ s/\W//g;
$rlon =~ s/\W//g;
$r->print('parameters missing'."\n"), return
unless $model and $variable and $rlat and $rlon;
my $sql = 'select name,number from dsm_season order by number';
my $sth = sql($r,$dbh,$sql);
my %seasons;
if ($sth) { for (1..$sth->rows) {
my ($name,$number) = $sth->fetchrow_array;
$seasons{$number-1} = $name;
}}
$r->print("$variable");
my $x = get_timeseries($r,$dbh,'seasonal',$model,$variable,'C','annual',$rlat,$rlon);
xml_scenario_timeseries($r,$x,'C',\%seasons);
if ($scenario) {
$x = get_timeseries($r,$dbh,'seasonal',$model,$variable,$scenario,'annual',$rlat,$rlon);
xml_scenario_timeseries($r,$x,$scenario,\%seasons);
}
$r->print("");
}
sub xml_scenario_timeseries {
my($r,$x,$scenario,$seasons) = @_;
unless (ref($x)) {
$r->print("$x: $scenario");
return;
}
$r->print("$scenario");
my %min;
my %max;
for (my $i = 0; $i < @$x; $i++) {
my $season = $seasons->{$i%4};
$min{$season} = $x->[$i] if !defined $min{$season} or $min{$season} > $x->[$i];
$max{$season} = $x->[$i] if !defined $max{$season} or $max{$season} < $x->[$i];
}
for my $season (keys %min) {
$r->print("$season$min{$season}$max{$season}");
}
$r->print("");
}
sub Distribution {
my($r,$q,$dbh,$u) = @_;
my $network = $q->param('network');
my $model = $q->param('model');
my $variable = $q->param('variable');
my $scenario = $q->param('scenario');
my $rlat = $q->param('rlat');
my $rlon = $q->param('rlon');
$network =~ s/[^\w\.-]//g;
$model =~ s/[^\w-]//g;
$variable =~ s/[^\w\s-]//g;
$scenario =~ s/[^\w-]//g;
$rlat =~ s/\W//g;
$rlon =~ s/\W//g;
$r->print('parameters missing'."\n"), return
unless $network and $model and $variable and $rlat and $rlon;
my @outcomes;
my $parser = new XML::DOM::Parser;
my $doc;
eval {
$doc = $parser->parsefile($conf->{dir}.'/dss/networks/'.$network);
};
if ($@ or !$doc) {
$r->print("There was an error in parsing $network: $@");
return;
} else {
my %variables;
my @variables = $doc->getElementsByTagName('VARIABLE');
for my $node (@variables) {
my @name_node = $node->getElementsByTagName('NAME');
next unless @name_node;
my $name = $name_node[0]->getFirstChild()->getNodeValue;
my $test_name = $name;
$test_name =~ s/[^\w\s-]//g;
next unless $test_name eq $variable;
my @list = $node->getElementsByTagName('OUTCOME');
for my $o (@list) {
push @outcomes, $o->getFirstChild()->getNodeValue;
}
}
}
$r->print("Variable $variable was not found in $network or it does not have outcome nodes.\n"), return
unless @outcomes;
my $sql = 'select name,number from dsm_season order by number';
my $sth = sql($r,$dbh,$sql);
my %seasons;
if ($sth) { for (1..$sth->rows) {
my ($name,$number) = $sth->fetchrow_array;
$seasons{$name} = $number;
}}
$sql = 'select v.name,v.acronym,u.name from dsm_variable v, dsm_unit u where v.unit=u.id';
$sth = sql($r,$dbh,$sql);
my %variables;
my %units;
if ($sth) { for (1..$sth->rows) {
my ($name,$acronym,$unit) = $sth->fetchrow_array;
$variables{$name} = $acronym;
$units{$acronym} = $unit;
}}
my ($var,$season,$comp_model) = parse_variable($variable,\%seasons,\%variables);
$r->print('could not parse variable$variable'."\n"), return
unless $var and $season;
my $unit = $units{$var};
my $percentage = 0;
my @bins;
for my $o (@outcomes) {
$o =~ s/\.\.\./ /;
my @n = $o =~ /([+-]?\d+\.\d+|[+-]?\d+\.|[+-]?\.\d+|[+-]?\d+)/g;
$r->print("Cannot parse outcome $o.\n"), return
unless @n == 2;
if ($o =~ /\%/) {
$percentage = 1;
for (@n) {
$_ /= 100;
}
} elsif ($unit eq 'mm/season' and $o =~ /cm\/season/) {
for (@n) {
$_ *= 10;
}
}
push @bins,@n;
}
my $type = $season eq 'annual' ? 'annual' : 'seasonal';
my $x;
if ($comp_model =~ /^deviation/) {
my $ctl = get_timeseries($r,$dbh,$type,$model,$var,'C',$season,$rlat,$rlon);
$r->print("no data\n"), return unless ref $ctl and @$ctl;
$x = get_timeseries($r,$dbh,$type,$model,$var,$scenario,$season,$rlat,$rlon);
$r->print("no data\n"), return unless ref $x;
my $average;
for (@$ctl) {
$average += $_;
}
$average /= @$ctl;
$r->print("couldn't compute deviation since average is zero\n"), return
unless $average;
for my $i (0..$#$x) {
if ($percentage) {
$x->[$i] = ($x->[$i]/$average)-1;
} else {
$x->[$i] = $x->[$i] - $average;
}
}
} else {
$x = get_timeseries($r,$dbh,$type,$model,$var,$scenario,$season,$rlat,$rlon);
$r->print("no data in $x\n"), return unless ref $x;
}
print STDERR "@bins\n@$x\n";
my $n = comp_distribution($x,@bins);
$r->print('');
$r->print("$variable");
$r->print("$network");
$r->print("$model");
$r->print("$var");
$r->print("$unit");
$r->print("$season");
$r->print("$comp_model");
$r->print("$scenario");
$r->print("$rlat");
$r->print("$rlon");
$r->print("",int($#bins/2)+1,"");
my $index = 0;
for (my $i = 0; $i < $#bins; $i += 2) {
$r->print("$index$bins[$i]",$bins[$i+1],"$n->[$i]");
$index++;
}
$r->print('');
}
sub parse_variable {
my($variable, $seasons, $variables) = @_;
my $var = '';
for my $name (keys %$variables) {
$var = $variables->{$name} if $variable =~ /$variables->{$name}/i;
}
unless ($var) {
for my $name (keys %$variables) {
$var = $variables->{$name} if $variable =~ /$name/i;
}
}
my $season = '';
for my $name (keys %$seasons) {
$season = $seasons->{$name} if $variable =~ /$name/i;
}
$season = 'annual' unless $season;
my $comp_model = '';
$comp_model = 'deviation from average in control' if $variable =~ /deviation/i;
return ($var,$season,$comp_model);
}
sub comp_distribution {
my ($X, @bins) = @_; # bins = min,max,min,max
my @n;
for my $x (@$X) {
for (my $i = 0; $i < $#bins; $i += 2) {
$n[$i]++ if $x >= $bins[$i] and $x < $bins[$i+1];
}
}
for (@n) {
$_ = 0 unless defined $_;
$_ /= @$X;
}
return \@n;
}
sub get_timeseries {
my ($r,$dbh,$type,$model,$variable,$scenario,$season,$rlat,$rlon) = @_;
if ($variable =~ /^frost/) {
my $sql = "select v from climate where var='$variable' and model='$model' and sc='$scenario' and rlat='$rlat' and rlon='$rlon'";
my $sth = sql($r,$dbh,$sql);
my @v;
if ($sth) { for (1..$sth->rows) {
my ($v) = $sth->fetchrow_array;
push @v,$v;
}}
return \@v;
}
my $fn = $conf->{dir}.'/dss/climate_data/'.$type.'.'.file($model,$variable,$scenario).'.nc';
my $ncobj;
eval {
$ncobj = PDL::NetCDF->new ($fn);
};
return $fn unless $ncobj;
my $varlist = $ncobj->getdimensionnames($variable);
# mokasin 1. kerralla netcdf:ien teon
$varlist = $ncobj->getdimensionnames('precip') unless @$varlist;
my @p1;
my @p2;
for my $v (@$varlist) {
if ($v eq 'time') {
push @p1,0;
push @p2,$ncobj->dimsize($v);
} elsif ($v eq 'rlat') {
push @p1,$rlat;
push @p2,1;
} elsif ($v eq 'rlon') {
push @p1,$rlon;
push @p2,1;
} elsif ($v eq 'height') {
push @p1,0;
push @p2,1;
}
}
my $slice;
eval {
$slice = $ncobj->get ($variable, \@p1, \@p2);
};
if ($@) {
$@ = '';
eval {
$slice = $ncobj->get ('precip', \@p1, \@p2);
};
}
return "$fn: $variable" if ($@);
my @x = list $slice;
return \@x if $season eq 'annual';
my @y;
for (my $i = $season-1; $i <= $#x; $i += 4) {
push @y,$x[$i];
}
return \@y;
}
sub file {
my ($model,$variable,$scenario) = @_;
my $d = dir($model,$variable,$scenario);
if ($model=~/^RCAO/) {
return "$variable.SMHI.$d";
} else {
my $y = $scenario eq 'C' ? '1960-1990' : '2070-2100';
return "$variable.HC.$d.$y";
}
}
sub dir {
my ($model,$variable,$scenario) = @_;
return 'HCCTL' if $model =~ /^RCAO-H/ and $scenario =~ /^C/;
return 'HCA2' if $model =~ /^RCAO-H/ and $scenario =~ /^A2/;
return 'HCB2' if $model =~ /^RCAO-H/ and $scenario =~ /^B2/;
return 'MPICTL' if $model =~ /^RCAO-E/ and $scenario =~ /^C/;
return 'MPIA2' if $model =~ /^RCAO-E/ and $scenario =~ /^A2/;
return 'MPIB2' if $model =~ /^RCAO-E/ and $scenario =~ /^B2/;
return 'adeha' if $model =~ /^Had/ and $scenario =~ /^C/;
return 'adhfa' if $model =~ /^Had/ and $scenario =~ /^A2/;
return 'adhfd' if $model =~ /^Had/ and $scenario =~ /^B2/;
die "$model,$variable,$scenario";
}