#!/usr/bin/perl -w

#This is the source code for the CLIME-DSS server.
#
#For more information see: clime.tkk.fi and
#http://geoinformatics.tkk.fi/twiki/bin/view/Main/CLIMEDSS
#
#The program code has been written by
#
#Ari Jolma ari.jolma at tkk.fi
#
#The program code is copyright Helsinki University of Technology
#
#This program is free software; you can redistribute it and/or modify
#it under the same terms as Perl itself. That is, using either
#
#	a) the GNU General Public License as published by the Free
#	Software Foundation; or
#
#	b) the "Artistic License".
#
#This program is distributed in the hope that it will be useful, but
#WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See either the
#GNU General Public License or the Artistic License for more details.
#
#You should have received a copy of the Artistic License with this Kit,
#in the file named "Artistic".  If not, get it from
#http://www.perl.com/language/misc/Artistic.html.
#
#You should also have received a copy of the GNU General Public License
#along with this program; if not, write to the Free Software
#Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
#USA.

BEGIN {
    $ENV{'PATH'} = '/bin:/usr/bin';
    delete @ENV{'IFS', 'CDPATH', 'ENV', 'BASH_ENV'};
}

use strict;

use CGI;
use DBI;
use UserAccount;
use FileHandle;
use File::stat;
use Date::Calc qw (:all);
use XML::DOM;
use ogr;
use Geo::Proj4;
use Geo::Raster;
use Geo::Vector;
use Gtk2::Ex::Geo::Renderer;
use PDL;
use PDL::NetCDF;

use vars qw($me $conf $ua $br $DENY $mime_types);

$me = 'wms.pl';
$DENY = 'Sorry, you are not authorized.';

use CLIME;
    
$conf = conf();

page();

sub page {

    $br = "<br>\n";

    my $r = Apache->request();
    my $q = new CGI;
    if ($conf->{error}) {
	serve_document($r,$q,{html=>'text/html'},$conf->{dir}.'/'.$conf->{error});
	return;
    }
    
    $ua = new UserAccount;
    
    my $x_cookie = $q->param('cookie');
    my $cookie = clime_ua($me,$r,$q,$ua,$conf,$x_cookie);

    $mime_types = mime_types($mime_types);

    my $service = my_param($q, 'SERVICE');
    my $request = my_param($q, 'REQUEST');
    my $doc = my_param($q, 'doc');
    my $cmd = my_param($q, 'cmd');

    if ($service eq 'file') {

	my $dir = "";
	if ($request =~ /^([\w\d\.]+)$/) {
	    $dir = $1;
	}
	my $file = "";
	if ($doc =~ /^([\w\d\.\-]+)$/) {
	    $file = $conf->{dir}."/dss/$dir/$1";
	}
	serve_document($r,$q,$mime_types,$file) if $file;

	return;

    }
    
#    unless (($x_cookie eq 'miracle') or ($ua->{currentuser} and $ua->{dbh})) {
    if (0) {

	$r->print($q->header( -type => 'text/html',
			      -cookie => $cookie ),
		  $q->start_html('Error'),
		  "You're not logged in or I've got no database.",
		  $q->end_html);

    } elsif ($service eq 'CLIME_DSS') {

	xml_header($r,$q);
	CLIME_DSS($r,$q,$ua);

    } elsif ($service eq 'WMS' and $request eq 'GetCapabilities') {

	xml_header($r,$q);
	GetCapabilities_1_1_1($r,$q);

    } elsif ($request eq 'GetMap') {

	GetMap_1_1_1($r,$q);
	
    } elsif ($_ = $doc or $cmd =~ /^Download/) {

	$_ = '' unless $_;

	xml_header($r,$q);
	
      SWITCH: {
	  MapSet($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^MapSet/;
	  LakeSet($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^LakeSet/;
	  LakeData($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^LakeData/;
	  ClimateVarsSet($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^ClimateVarsSet/;
	  ClimateTimeseries($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^ClimateTimeseries/;
	  TopicSet($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^TopicSet/;
	  Arrow($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^Arrow/;
	  Climate($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^Climate/;
	  Patterns($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^Patterns/;
	  Distribution($r,$q,$ua->{dbh},$ua->{currentuser}), last SWITCH  if /^Distribution/;
      }
	
	my $cmd = $q->param('cmd') || '';
	if ($cmd =~ /^Download/) {
	    my $network = $q->param('network');
	    my $filename = $network;
	    if (open(XML,$conf->{dir}.'/dss/networks/'.$filename)) {
		while (<XML>) {
		    $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('<?xml version="1.0" encoding="ISO-8859-1"?>'."\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('<p>',$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,'</p>');
}

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.<br>\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("<p>Now documenting network '$network' (",href('back',$me),")</p>\n",'<p>',
		  $q->start_form,
		  $q->hidden('network',$network));

	my $descr = exists $mapping->{descr} ? $mapping->{descr} : '';

	$r->print("Description of <b>$network</b>:<br>\n",
		  $q->textarea(-name=>"descr",
			       -default=>$descr,
			       -rows=>5,
			       -columns=>40),"<br>\n");
	
	my $topic = exists $mapping->{topic} ? $mapping->{topic} : 'Nothing';
	
	$r->print("The <b>topic</b> of this network is ",
		  $q->popup_menu(-name=>"topic",
				 -values=>\@topics,
				 -default=>$topic),"<br><br>\n");
	
	my $csite = exists $mapping->{csite} ? $mapping->{csite} : 'Nothing';
	
	$r->print("The <b>calibration site</b> of this network is ",
		  $q->popup_menu(-name=>"csite",
				 -values=>\@csites,
				 -default=>$csite),"<br><br>\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 <b>$variable</b>:<br>\n",
		      $q->textarea(-name=>"descr var $i",
				   -default=>$descr,
				   -rows=>5,
				   -columns=>40),"<br>\n");


	    $r->print("Link <b>$variable</b> to ",
		      $q->popup_menu(-name=>"var $i",
				     -values=>$distributions,
				     -default=>$default),"<br>\n");

	}
	
	$r->print('<p>',$q->submit('cmd','Store this'),'</p>',$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:<br/><font color=\"red\">$@</font><br/>");
	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("<TopicSet>");

    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("<Topic TopicCode=\"$TopicCode\">");
	my $n = $network;
	$n =~ s/\.xml$//i;
	my @n = split /_/,$n;
	if ($mappings->{$network}{topic}) {
	    $r->print("<Name>$mappings->{$network}{topic}</Name>");
	} else {
	    $r->print("<Name>$n[0]</Name>");
	}
	$r->print("<Site>$n[1]</Site>");
	$r->print("<Description>$mappings->{$network}{descr}</Description>") if $mappings->{$network}{descr};
	$r->print("<NetworkXML>$network</NetworkXML>");
	for my $model ('RCAO-E','RCAO-H','HadRM3p') {
	    my $g = $model;
	    $g =~ s/\-.*$//;
	    my $lake = $mappings->{$network}{csite};
	    $r->print("<NoCalibrationSite model=\"$model\"/>") 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("<SuitabilityMap model=\"$model\">",
			  "$conf->{URLRoot}$me?SERVICE=file\&amp\;REQUEST=sims\&amp\;doc=$fn",
			  "</SuitabilityMap>");
	    } else {
		$r->print("<MissingSuitabilityMap model=\"$model\">$file</MissingSuitabilityMap>");
	    }
	}
#	$r->print("<SuitabilityMap model=\"RCAO-E\">$conf->{URLRoot}?dss=sims/RE_6065.png</SuitabilityMap>");
#	$r->print("<SuitabilityMap model=\"RCAO-H\">$conf->{URLRoot}?dss=sims/RH_6065.png</SuitabilityMap>");
#	$r->print("<BNURI>$conf->{URLRoot}?dss=networks/$network</BNURI>");
	$r->print("<BNURI>$conf->{URLRoot}$me?network=$network\&amp\;cmd=Download+the+model</BNURI>");

	my $parser = new XML::DOM::Parser;
	my $doc;
	eval {
	    $doc = $parser->parsefile($conf->{dir}.'/dss/networks/'.$network);
	};
	if ($@ or !$doc) {
	    $r->print("<error type=\"parsing\"><network>$network</network>$@</error>");
	    $r->print("</Topic>");
	    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("<Variable Name=\"$v\">");

		    my $vdata = $mappings->{$network}{$v};

		    my $descr = $vdata->{descr} if $vdata;
		    $r->print("<Description>$descr</Description>");

		    my $sp_dist = $vdata->{sp_dist} if $vdata;

		    unless ($sp_dist) {
			$r->print("<DataURL/></Variable>");
			next;
		    }

		    $r->print("<SpatialVariable>$sp_dist</SpatialVariable>");

		    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("<Spatial_variable>$sp_dist</Spatial_variable>");
		    $r->print("<prefix>$prefix</prefix>");

		    my $count = 0;

		    if (0) { # not used any more
		    for my $fn (@$distributions) {

#			$r->print("<fn>$fn</fn>");

			for my $model (@ClimateModels)  {
			    for my $scenario (@ClimateScenarios) {

#				if ($fn =~ /$sp_dist/) {
				if ($fn eq $prefix.$model.'_'.$scenario.'_'.$sp_dist.'.png') {

				    $r->print("<DataURL Scenario=\"$scenario\" Model=\"$model\">",
					      "$conf->{URLRoot}?dss=distributions/$fn</DataURL>");
				    $count++;
				    
				}
			    }
			}
		    }
		    $r->print("<DataURL/>") unless $count;
		}
		    $r->print(" </Variable>");
		}
	    }
	}

	$doc->dispose;

	$r->print("</Topic>");
    }
    $r->print("</TopicSet>");
}

sub ClimateVarsSet {
    my($r,$q,$dbh,$u) = @_;
    $r->print("<ClimateVarsSet>");

    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("<ClimateScenarioSet>");
    my $id = 1;
    for my $scenario (@ClimateScenarios) {
#	next if $scenario eq 'C'; # needed!
	$r->print("<ClimateScenario ID=\"$id\">");
	$ClimateScenarioID{$scenario} = $id;
	$id++;
	$r->print("<Name>$scenario</Name>");
	$r->print("<Description>$scenarioHelp{$scenario}</Description>");	    	    
	$r->print("</ClimateScenario>");
    }
    $r->print("</ClimateScenarioSet>");
    
    my %ClimateModelID;
    $r->print("<ClimateModelSet>");
    $id = 1;
    for my $model (@ClimateModels) {
    	    $r->print("<ClimateModel ID=\"$id\">");
    	    $ClimateModelID{$model} = $id;
    	    $id++;
	    $r->print("<Name>$model</Name>");
	    my $grid = $model =~ /^Had/ ? 'HadRM3p_s' : 'RCAO';
    	    $r->print("<GridURL>$conf->{URLRoot}$me?SERVICE=file\&amp\;REQUEST=grids\&amp\;doc=${grid}_g.png</GridURL>");
    	    $r->print("<Description>$modelsHelp{$model}</Description>");
    	    $r->print("</ClimateModel>");
    }
    $r->print("</ClimateModelSet>\n");

    my %SeasonID;
    $r->print("<SeasonSet>");
    $id = 1;
    for my $season (@seasons) {
	$r->print("<Season ID=\"$id\">");
	$SeasonID{$season} = $id;
	$id++;
	$r->print("<Name>$season</Name>");
	$r->print("<Description>$seasonHelp{$season}</Description>");
	$r->print("</Season>");
    }
    $r->print("</SeasonSet>\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("<MapTypeSet>");
        $id = 1;
        for my $mapType (@MapTypes) {
		$r->print("<MapType ID=\"$id\">");
		$MapTypeID{$mapType} = $id;
		$id++;
		$r->print("<Name>$mapType</Name>");
		$r->print("<Description>$mapTypeHelp{$mapType}</Description>");
		$r->print("</MapType>");
        }
    $r->print("</MapTypeSet>\n");

    my %VariableID;
    $r->print("<VariableSet>");
    $id = 1;
    for my $variable (@variables) {
	$r->print("<Variable ID=\"$id\">");
	$VariableID{$variable} = $id;
	$id++;
	$r->print("<Acronym>$variable</Acronym>\n");
	$r->print("<Scale>$scale{$variable}</Scale>\n");
	$r->print("<Unit>$unit{$variable}</Unit>\n");
	$r->print("<Name>$fname{$variable}</Name>\n");
	$r->print("<Description>$variableHelp{$variable}</Description>\n");

	for my $rule (@{$rules{$variable}}) {
	    for my $key (keys %$rule) {
		my $k = $key;
		$k =~ s/\</&lt;/g;
		$k =~ s/\>/&gt;/g;
		$r->print("<Rule category=\"$rule->{$key}\">$k</Rule>\n");
	    }
	}
	
	$r->print("</Variable>\n");
    }
    $r->print("</VariableSet>\n");

    my $URLRoot = $conf->{URLRoot};

    for my $variable (@variables) { 
	$r->print('<ClimateVariable VariableID="',$VariableID{$variable},'">');	
	for my $scenario (@ClimateScenarios) { 	
		for my $model (@ClimateModels) {

		    for my $season (@seasons) {
			my $fn = $model.'_'.$scenario.'_'.lc($season).'_'.$variable.'.'.$ext{1};
			$r->print("<VisualizationURL SeasonID=\"$SeasonID{$season}\" ClimateScenarioID=\"$ClimateScenarioID{$scenario}\" ClimateModelID=\"$ClimateModelID{$model}\"  MapTypeID=\"1\">",
				  $URLRoot,$me,"?SERVICE=file\&amp\;REQUEST=",$subdir{2},"\&amp\;doc=",$fn,
				  "</VisualizationURL>");
		    }
		    
		    for my $season (@seasons) {
			my $fn = $model.'_'.$scenario.'_'.lc($season).'_'.$variable.'.'.$ext{1};
			$r->print("<VisualizationURL SeasonID=\"$SeasonID{$season}\" ClimateScenarioID=\"$ClimateScenarioID{$scenario}\" ClimateModelID=\"$ClimateModelID{$model}\"  MapTypeID=\"2\">",
				  $URLRoot,$me,"?SERVICE=file\&amp\;REQUEST=",$subdir{3},"\&amp\;doc=",$fn,
				  "</VisualizationURL>");
		    }		    
		}
	}
	$r->print("</ClimateVariable>");
    }
    
    $r->print("</ClimateVarsSet>");
}

sub MapSet {
    my($r,$q,$dbh,$u) = @_;
    $r->print('<MapSet>'."\n");

    $r->print('<Map>'."\n");
    $r->print("<ID>1</ID>");
    $r->print("<Name>Physical (relief)</Name>");
    $r->print("<WMS>No</WMS>");
    $r->print("<URL>$conf->{URLRoot}?dss=physical_map.jpg</URL>");
    $r->print("<Description>On this relief map the colors are in following order ",
	      "from lowest to highest: blue, green, yellow, brown, purple.</Description>");
    $r->print('</Map>'."\n");

    $r->print('<Map>'."\n");
    $r->print("<ID>2</ID>");
    $r->print("<Name>Land cover (CORINE)</Name>");
    $r->print("<WMS>Yes</WMS>");
    my $m = $me;
    $m =~ s/^\///;
    $r->print("<URL>$conf->{URLRoot}$m?map=corine</URL>");
    $r->print("<Description>Land cover map is based on CORINE data.</Description>");
    $r->print('</Map>'."\n");

    $r->print('<Map>'."\n");
    $r->print("<ID>3</ID>");
    $r->print("<Name>Peatlands and lakes</Name>");
    $r->print("<WMS>No</WMS>");
    $r->print("<URL>$conf->{URLRoot}?dss=peat_and_lakes.png</URL>");
    $r->print("<Description>The peatlands are marked with brown color. ",
	      "The darker the colour, the more peat there is. ",
	      "Same peat land data is used in DOC modeling. </Description>");
    $r->print('</Map>'."\n");
    
    $r->print('</MapSet>'."\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('<LakeSet>'."\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("<Lake SiteCode=\"$sitecode\">\n");
	$r->print("<DBID>$id</DBID>\n");
	$r->print("<Name>$lake</Name>\n");
#	$r->print("<Location><X>$ij[1]</X><Y>$ij[0]</Y></Location>\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('</Lake>'."\n");
	$sitecode++;

    }}

    $r->print('</LakeSet>'."\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('<Error>parameters missing</Error>'."\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('<Arrow>');
    $r->print("<sql>$sql</sql>");
    $r->print("<Model>$model</Model>");
    $r->print("<Scenario>$scenario</Scenario>");
    $r->print("<from_rlat>$rlat</from_rlat>");
    $r->print("<from_rlon>$rlon</from_rlon>");
    $r->print("<to_rlat>$from_rlat</to_rlat>");
    $r->print("<to_rlon>$from_rlon</to_rlon>");
    $r->print("<from_i>$from_i</from_i>");
    $r->print("<from_j>$from_j</from_j>");
    $r->print("<to_i>$to_i</to_i>");
    $r->print("<to_j>$to_j</to_j>");
    $r->print("<similarity>$sim</similarity>");
    $r->print('</Arrow>'."\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('<Error>parameters missing</Error>'."\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('<Climate>');
    $r->print("<Model> $model </Model>");
    $r->print("<Season> $season </Season>");
    $r->print("<rlat> $rlat </rlat>");
    $r->print("<rlon> $rlon </rlon>");
#    $r->print("<sql sql=\"$sql\"/>");
    my $sc0 = '';
    if ($sth) { 
	for (1..$sth->rows) {
	    my ($sc,$var,$v) = $sth->fetchrow_array;
	    if ($sc ne $sc0) {
		if ($sc0) {
		    $r->print("<VariableSeasonalAverage Name=\"snow\">no data</VariableSeasonalAverage>") if $model eq 'HadRM3p';
		    $r->print('</Scenario>');
		}
		$r->print("<Scenario Name=\"$sc\">");
		$sc0 = $sc;
	    }
	    $r->print("<VariableSeasonalAverage Name=\"$var\"> $v</VariableSeasonalAverage>");
	}
	$r->print("<VariableSeasonalAverage Name=\"snow\">no data</VariableSeasonalAverage>") if $model eq 'HadRM3p';
	$r->print('</Scenario>') if $sth->rows;
    }
    $r->print('</Climate>');
}

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('<Error>parameters missing</Error>'."\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('<Patterns>');
    #$r->print("<sql> $sql </sql>");
    $r->print("<Model> $model </Model>");
    $r->print("<Season> $season </Season>");
    $r->print("<rlat> $rlat </rlat>");
    $r->print("<rlon> $rlon </rlon>");

    my @scenarios = ('C','A2','B2');

    my $i = 1;
    for my $pattern (@patterns) {
	$r->print("<Pattern ID=\"$i\">");
	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("</Pattern>");
    }

    $r->print('</Patterns>'."\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('<Error>parameters missing</Error>'."\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("<ClimateTimeseries><Variable>$variable</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("</ClimateTimeseries>");
    
}

sub xml_scenario_timeseries {
    my($r,$x,$scenario,$seasons) = @_;
    unless (ref($x)) {
	$r->print("<Error>$x: $scenario</Error>");
	return;
    }

    $r->print("<Scenario><Name>$scenario</Name>");
    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("<Value><Season>$season</Season><Min>$min{$season}</Min><Max>$max{$season}</Max></Value>");
    }
    $r->print("</Scenario>");
}

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('<Error>parameters missing</Error>'."\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("<Error>There was an error in parsing $network: $@</Error>");
	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("<Error>Variable $variable was not found in $network or it does not have outcome nodes.</Error>\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('<Error>could not parse variable<variable>$variable</variable></Error>'."\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("<Error>Cannot parse outcome $o.</Error>\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("<Error>no data</Error>\n"), return unless ref $ctl and @$ctl;
	$x = get_timeseries($r,$dbh,$type,$model,$var,$scenario,$season,$rlat,$rlon);
	$r->print("<Error>no data</Error>\n"), return unless ref $x;
	my $average;
	for (@$ctl) {
	    $average += $_;
	}
	$average /= @$ctl;
	$r->print("<Error>couldn't compute deviation since average is zero</Error>\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("<Error>no data in $x</Error>\n"), return unless ref $x;
    }
    print STDERR "@bins\n@$x\n";

    my $n = comp_distribution($x,@bins);
    
    $r->print('<Distribution>');

    $r->print("<NAME>$variable</NAME>");
    $r->print("<network>$network</network>");
    $r->print("<model>$model</model>");
    $r->print("<variable>$var</variable>");
    $r->print("<unit>$unit</unit>");
    $r->print("<season>$season</season>");
    $r->print("<comp_model>$comp_model</comp_model>");
    $r->print("<scenario>$scenario</scenario>");
    $r->print("<rlat>$rlat</rlat>");
    $r->print("<rlon>$rlon</rlon>");

    $r->print("<bins>",int($#bins/2)+1,"</bins>");
    my $index = 0;
    for (my $i = 0; $i < $#bins; $i += 2) {
	$r->print("<bin><index>$index</index><from>$bins[$i]</from><to>",$bins[$i+1],"</to><probability>$n->[$i]</probability></bin>");
	$index++;
    }

    $r->print('</Distribution>');
}

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";
}