package EnhancedGrid;
use strict;
use warnings;
use SimpleObject;
use Storable qw(dclone);
use CtUtilities qw(round_to is_in_list);
use MasRaster;
use Carp;
our @ISA = qw(SimpleObject);


sub create {
	my $class = shift;	
	croak "Error: create is a class method of class EnhancedGrid and cannot be called as an object method!\n" if (ref $class);
	my @attributes = @_;
	my %precheck;
	while (@attributes) {
		my $r = shift (@attributes);
		my $s = shift (@attributes);
    		$precheck{$r} = $s;
	}
	#Defaults
  	$precheck{ncols} = 1 if (not exists $precheck{ncols});
  	$precheck{nrows} = 1 if (not exists $precheck{nrows});
  	$precheck{NoData} = -1 if (not exists $precheck{NoData});
	
  	$precheck{xllcorner} = 0 if (not exists $precheck{xllcorner});
	$precheck{yllcorner} = 0 if (not exists $precheck{yllcorner});
	$precheck{cellsize} = 1 if (not exists $precheck{cellsize});
	my $cols = $precheck{ncols};
	my $rows = $precheck{nrows};
	my $fd = (exists $precheck{FirstData}) ? $precheck{FirstData} : $precheck{NoData};
	delete $precheck{FirstData} if (exists $precheck{FirstData});
	$precheck{Raster} = MasRaster->create($cols, $rows, $fd);
	my $self = $class->SUPER::create(%precheck);
	
	return $self;
}

sub create_alike {
	my $self = shift;
	my $filled_with = shift;
	my $class = ref($self);
	croak "You cannot use create_alike as a class method! (called from $self)" unless (ref $self);
	my %precheck;
	for my $i (qw/ncols nrows xllcorner yllcorner cellsize NoData/) {
		$precheck{$i} = $self->get_attribute($i);
	}
	$precheck{FirstData} = (defined $filled_with) ? $filled_with: $precheck{NoData};
	
	my $newinstance = $class->create(%precheck);
	
}
sub create_map_from_sequence {
	my ($class, $arrayref, $mapheight, $nodata, $name, $cellsize) = @_;
	my @array = @$arrayref;
	#append nodata cells to fill otherwise missing cells
	push @array, ($nodata) x ($mapheight - (scalar(@array) % $mapheight)) if (scalar(@array) % $mapheight > 0) ;
	
	my @map;
	for my $i (0 .. scalar (@array) / $mapheight - 1) {
	  push @map, [ @array[$i * $mapheight .. ($i+1) * $mapheight -1] ];
	}
	 my %precheck;
	$precheck{ncols} = scalar(@map) ;
  	$precheck{nrows} = $mapheight;
  	$precheck{NoData} = $nodata;
	
  	$precheck{xllcorner} = 0 ;
	$precheck{yllcorner} = 0;
	$precheck{cellsize} = $cellsize;
	$precheck{id} = $name;	
	$precheck{Raster} = MasRaster->create_from_preinitialized(\@map);
	my $self = $class->SUPER::create(%precheck);
	return $self;
}

sub ncols {	
	my $self = shift;
 	return $self->get_attribute("ncols");
}
sub nrows {	
	my $self = shift;
 	return $self->get_attribute("nrows");
}

sub move_to {
	my $self = shift;	
	my $x = shift;
	my $y = shift;
	$self->set_attribute("xllcorner", $x);
	$self->set_attribute("yllcorner", $y);	
}
sub get_cell {
	my $self = shift;
	my $x = shift;
	my $y = shift;
	$self->{Raster}->get_cell($x,$y);
}
sub set_cell {
	my $self = shift;
	my $x = shift;
	my $y = shift;
	my $value = shift;
	$self->{Raster}->set_cell($x,$y, $value);
}
sub create_polygonraster {
	my $self = shift;
	my $onlythesevalues = shift;
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	$self->{PolygonRaster} = MasRaster->create($ncols, $nrows, -1);
	$self->{PolygonTable} = [];
	
	$self->refresh_polygonraster($onlythesevalues);
}

sub refresh_polygonraster {
	my $self = shift;
	my $onlythesevalues = shift;
	#check whether  Polygon raster exists
	croak "Error: Cannot refresh PolygonRaster that has not been created yet!\n" unless (defined $self->{PolygonRaster});
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	$self->{PolygonRaster}->fill(-1);
	$self->{PolygonTable} = undef;
	my $curPolygon = -1;
	# loop over all cells of the grid for possible start cells
	for my $startX (0..$ncols -1) {
		for my $startY (0..$nrows -1) {
			#do not start scan if cell is already part of a polygon
			next if ($self->{PolygonRaster}->get_cell($startX, $startY) > -1);

			#do not start if polygons of these values shall not be included
			my $v = $self->{Raster}->get_cell($startX, $startY);
			if (defined $onlythesevalues) {
				if (ref($onlythesevalues) eq "ARRAY") {
					next unless is_in_list($v, $onlythesevalues );
				}	
				elsif (not ref $onlythesevalues ) {
					next unless $v == $onlythesevalues;
				}
			}
			#initialize new polygon
			$self->new_polygon($startX, $startY,undef)
		}
	}
	return 1;
}
sub new_polygon {
	my $self = shift;
	my $startX = shift;
	my $startY = shift;
	my $polygon_id = shift;
	$polygon_id = $#{ $self->{PolygonTable} } + 1 unless (defined $polygon_id);
	#$polygon_id = 0 if (not defined ($polygon_id) or $polygon_id < 0 );
	#print "Create polygon $polygon_id\n";
	$self->{PolygonTable}->[$polygon_id]->{value} = $self->{Raster}->get_cell($startX, $startY);
	$self->{PolygonTable}->[$polygon_id]->{handlex} = $startX;
	$self->{PolygonTable}->[$polygon_id]->{handley} = $startY;	
	$self->find_and_attach_cells_to_polygon($polygon_id);
	return $polygon_id;
}
sub delete_polygon {
        my $self = shift;
	my $p = shift;
	$self->{PolygonTable}->[$p]->{size} = -1;
        $self->{PolygonTable}->[$p]->{value} = -9999;
	$self->{PolygonTable}->[$p]->{handlex} = -1;
	$self->{PolygonTable}->[$p]->{handley} = -1;
}
#
sub find_and_attach_cells_to_polygon {
	my $self = shift;
	my $polygon_id = shift;
	$self->{PolygonTable}->[$polygon_id]->{size} = 0;
	#start scanning cells with start cells, 
	#iterative because recursive does reach limits in perl
	my $ubound = $self->get_attribute("nrows");
	my $rbound = $self->get_attribute("ncols");
	my @stack;
	push @stack, {x =>$self->{PolygonTable}->[$polygon_id]->{handlex}, y => $self->{PolygonTable}->[$polygon_id]->{handley}, from => "start"};
	while (@stack) {
		my $c = pop @stack;
		my $curX = $c->{x};
		my $curY = $c->{y};
		my $from = $c->{from};
		
		if ($self->{Raster}->get_cell($curX,$curY) ==  $self->{PolygonTable}->[$polygon_id]->{value} and $self->{PolygonRaster}->get_cell($curX,$curY) != $polygon_id) {
			$self->{PolygonRaster}->set_cell($curX,$curY, $polygon_id);
			$self->{PolygonTable}->[$polygon_id]->{size} += 1;
			push @stack, {x => $curX - 1, y=> $curY, from => "right"}   if ($curX -1 >= 0 and not $from eq "left");
			push @stack, {x => $curX + 1, y=> $curY, from => "left"}   if ($curX +1 < $rbound and not $from eq "right");
			push @stack, {x => $curX, y=> $curY - 1, from => "above"}   if ($curY -1 >= 0 and not $from eq "below");
			push @stack, {x => $curX, y=> $curY + 1, from => "below"}   if ($curY +1 < $ubound and not $from eq "above");
		}
		
	}
	return 1;
}
sub update_polygonraster {
	my $self = shift;
	my @cells_changed = @{ shift @_ };
	my $new_value = shift;
	my $level = shift; # 0 -> only reduce size of existing polygons, create new polygons only if separated from old polygon, 1 -> reduce size of existing polygons and create new polygons for changed cells, 2 -> adapt (reduce or enlarge) size of existing polygons and create new ones
				#only 0 implemented so far
	##make list of cell names
	my @cell_changed_reflist = map {$_->{x}."_".$_->{y}  } @cells_changed;
	
	##reduce size of existing polygons
	my %affected_polygons; #hash instead of array to ensure uniqueness
	my @neighboring_cells;
	for my $c (@cells_changed) {
		my $p = $self->{PolygonRaster}->get_cell($c->{x},$c->{y});
		if ($self->{PolygonTable}->[$p]->{value} != $new_value )  {
			if (exists $affected_polygons{$p} ){
				$affected_polygons{$p}->{cells_lost} += 1 ;
			}
			else {
				$affected_polygons{$p}->{cells_lost} = 1 ;
				$affected_polygons{$p}->{neighboring_cells_to_test} = {};
			}
			$affected_polygons{$p}->{handle_cell_affected} = 1 if ($c->{x} == $self->{PolygonTable}->[$p]->{handlex} and $c->{y} == $self->{PolygonTable}->[$p]->{handley});
			for my $d ( $self->neighborhood($c ,4) ) {
				my $cellref = $d->{x}."_".$d->{y};
				$affected_polygons{$p}->{neighboring_cells_to_test}->{$cellref} = 1  unless (is_in_list($cellref, \@cell_changed_reflist) 
					or $p != $self->{PolygonRaster}->get_cell($d->{x}, $d->{y}) )  ;	
			}
			#set polygon_id of cell to -1! -- can be updated later if new polygons are created
			$self->{PolygonRaster}->set_cell($c->{x},$c->{y},-1);
		}
	}
	for my $p (keys %affected_polygons) {
		if ($self->{PolygonTable}->[$p]->{size} == $affected_polygons{$p}->{cells_lost} ) {
			$self->delete_polygon($p);
		    	next;
		}
		my @cells = keys %{ $affected_polygons{$p}->{neighboring_cells_to_test} };
		my @start_cells_for_new_polygon;
		while (@cells) {
			my $c = shift @cells;
			$c =~ /(\d+)_(\d+)/;
			my $c_cell = {x => $1, y=> $2};
			push  @start_cells_for_new_polygon, $c_cell;
			#print "\n";
			my @cells_temp;
			while (@cells) {
				my $d = shift @cells;
				$d =~ /(\d+)_(\d+)/;
				my $d_cell = {x => $1, y=> $2};
				unless ($self->a_connection_exists_between($c_cell, $d_cell)) {
					push @cells_temp, $d;
				}
			};
			@cells = @cells_temp;
			
		}
		#if @start_cells_for_new_polygon contains only one cell, the polygon has not been split
		if (scalar @start_cells_for_new_polygon == 1) {
			$self->{PolygonTable}->[$p]->{size} -= $affected_polygons{$p}->{cells_lost};
			if ($affected_polygons{$p}->{handle_cell_affected}) {
				$self->{PolygonTable}->[$p]->{handlex} = $start_cells_for_new_polygon[0]->{x};
				$self->{PolygonTable}->[$p]->{handley} = $start_cells_for_new_polygon[0]->{y};
			}
		}
		elsif (scalar @start_cells_for_new_polygon > 1) {
			#map {print "$_->{x}, $_->{y}\n"; } @start_cells_for_new_polygon;
			#first one becomes start cell of old polygon
			my $c = shift @start_cells_for_new_polygon; 
			$self->{PolygonTable}->[$p]->{handlex} = $c->{x};
			$self->{PolygonTable}->[$p]->{handley} = $c->{y};
			#print "Start cell old polygon $c->{x}, $c->{y}\n";
			#others become new polygons
			my $size_in_others = 0;
			foreach my $d (@start_cells_for_new_polygon) {
				#print "Creating new polygon with start cell $d->{x}, $d->{y}\n";
				my $newid= $self->new_polygon($d->{x}, $d->{y}, undef);
				$size_in_others += $self->{PolygonTable}->[$newid]->{size};
			}
			
			#adapt size  of "old" polygon
			$self->{PolygonTable}->[$p]->{size} -= $size_in_others + $affected_polygons{$p}->{cells_lost};
			
		}
	}
	
}



sub print_polygontable {
	my $self = shift;
	croak "Error: No Polygonraster available for this grid!\n" unless (exists $self->{PolygonRaster});
	for my $i (0..$#{$self->{PolygonTable}}) {
		print "$i\t$self->{PolygonTable}->[$i]->{value}\t$self->{PolygonTable}->[$i]->{size}\t$self->{PolygonTable}->[$i]->{handlex},$self->{PolygonTable}->[$i]->{handley}\n";
	}
}

sub get_id_and_size_of_largest_polygon_of_value {
	my $self	=	shift;
	my $valuesought	=	shift;
	my $id_largest  = undef;
	my $size_largest = -1;
	for my $i (0..$#{$self->{PolygonTable}} ) {
		if ($self->{PolygonTable}->[$i]->{value} == $valuesought and $self->{PolygonTable}->[$i]->{size} > $size_largest) {
			$id_largest = $i;
			$size_largest = $self->{PolygonTable}->[$i]->{size};
		}
	}
	return ($id_largest, $size_largest);
}
sub get_id_of_first_polygon_of_value_and_larger_or_equal {
	my $self	=	shift;
	my $valuesought	=	shift;
	my $minsize	= shift;
	my $id_first  = undef;
	for my $i (0..$#{$self->{PolygonTable}} ) {
		if ($self->{PolygonTable}->[$i]->{value} == $valuesought and $self->{PolygonTable}->[$i]->{size} >= $minsize) {
			$id_first = $i;
			last;
		}
	}
	return $id_first;
	
}

sub get_id_of_random_polygon_of_value_and_larger_or_equal {
	my $self	=	shift;
	my $valuesought	=	shift;
	my $minsize	=	shift;
	my @options;
	#print "Valuesought: $valuesought , Minsize: $minsize\n";
	for my $i (0..$#{$self->{PolygonTable}} ) {
		if ($self->{PolygonTable}->[$i]->{value} == $valuesought and $self->{PolygonTable}->[$i]->{size} >= $minsize) {
			push @options, $i;
		}
	}
	my $r = int rand(@options);
	return $options[$r];
}

sub get_handlecell_of_polygon_by_id {
	my $self = shift;
	my $id	= shift;
	my $x = $self->{PolygonTable}->[$id]->{handlex};
	my $y = $self->{PolygonTable}->[$id]->{handley};
	return ($x,$y);

}

sub get_polygon_id_by_cell {
	my $self = shift;
	my $x = shift;
	my $y = shift;
	return $self->{PolygonRaster}->get_cell($x,$y);
}
sub get_size_of_polygon_id {
	my $self = shift;
	my $id = shift;
	$self->{PolygonTable}->[$id]->{size}
}
sub get_value_of_polygon_id {
	my $self = shift;
	my $id = shift;
	$self->{PolygonTable}->[$id]->{value}
}
sub get_polygon_size_by_cell {
	my $self = shift;
	my $x = shift;
	my $y = shift;
	my $id = $self->{PolygonRaster}->get_cell($x,$y);
	$self->{PolygonTable}->[$id]->{size};	
}
sub neighborhood {
	my $self = shift;
	$self->{Raster}->neighborhood(@_);
}
sub a_connection_exists_between {
	my $self = shift;
	$self->{Raster}->a_connection_exists_between(@_);
}
sub print_valuetable {
	my $self = shift;
	croak "Error: No Valuetable available for this grid!\n" unless (exists $self->{ValueTable});
	my @Sorted = sort {$a <=> $b } keys %{$self->{ValueTable}};
	foreach (@Sorted) {
		print "$_\t$self->{ValueTable}->{$_}\n";
	}
	
}

sub update_valuetable {
	my $self = shift;
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my %valuetab;
	for my $x (0..$ncols -1) {
		for my $y (0..$nrows -1) {
			my $value = $self->get_cell($x,$y);
			$valuetab{$value} = 0 unless (exists $valuetab{$value});
			$valuetab{$value} += 1;
		 }
	}
	$self->{ValueTable} = \%valuetab;
}

sub save_polygontable_as_ascii {
	my $self = shift;
	my $fname = shift;
	croak "Error: No Polygonraster available for this grid!\n" unless (exists $self->{PolygonRaster});
	open (my $fh, ">", $fname) or croak "$!";
	print $fh "PolygonID\tValue\tSize\tX\tY\n";
	for my $i (0..$#{$self->{PolygonTable}}) {
		print $fh "$i\t$self->{PolygonTable}->[$i]->{value}\t$self->{PolygonTable}->[$i]->{size}\t$self->{PolygonTable}->[$i]->{handlex}\t$self->{PolygonTable}->[$i]->{handley}\n";
	}
	close $fh;
}


sub save_as_ascii {
	my $self = shift;
	my $fname = shift;
	my $what = shift;
	
	if (defined $what ) {
		my @listofdefinedrasters = ("Raster", "PolygonRaster");
		croak "Error: $what is not a grid layer of the EnhancedGrid Class" unless (grep $_ eq $what, @listofdefinedrasters);

	}
	else {
		$what = "Raster"
	}
	#print "Saving $what as $fname ...\n";
	open (my $fh, ">", $fname) or croak "Error when trying to open $fname: $!";
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my $xllcorner	= $self->get_attribute("xllcorner");
	my $yllcorner	= $self->get_attribute("yllcorner");
	my $cellsize	= $self->get_attribute("cellsize");
	my $NODATA_value	= $self->get_attribute("NoData");
	
	print $fh <<HD;
ncols		$ncols
nrows		$nrows
xllcorner	$xllcorner
yllcorner	$yllcorner
cellsize	$cellsize
NODATA_value	$NODATA_value
HD

	for my $j (reverse (0..$nrows-1)  ){
		for my $i (0..$ncols-1){
			
			if (defined $self->{$what}->get_cell($i,$j)) {
				print $fh $self->{$what}->get_cell($i,$j)." ";
			} 
			else {
				print "Warning: Empty cell: $i, $j\n";
			}
		}
		print $fh "\n";

	}
	
	close $fh;
}

sub create_from_ascii {
	my $class = shift;
	croak "Error: Create_from_ascii is a class method, you have called it from an object!\n" if (ref $class);
	my $file = shift;
	
	#print "Reading data from $file...\n";
	open (my $infile, "<", $file) or croak "Error when opening $file: $! \n\n";
	my $i = 0;
	my $currow = 0;
	my %attributes;
	my @PreRaster;
	while (<$infile>) {
		chomp;
		my $Line =  $_;
		next if ($Line =~ /^\s*$/); #skip empty lines
		$Line =~ s/^\s+//; # cut leading spaces to be robust against formating errors
		$i++;
		
		#read compulsory header
		if ($i <=5) {
			my ($key, $value) = split '\s+', $Line, 2;
			$attributes{lc($key)} = $value;
#			print "$key => $value\n";
		}
		#read nodata-value if present
		elsif ($i == 6 and lc($Line) =~ /^\s*"?nodata/) {
			my ($key, $value) = split /\s+/, $Line, 2;
			$attributes{"NoData"} = $value;
#			print "$key => $value\n";
		}
		#read cell values
		else {
			my @row = split /\s+/, $Line;
			push @PreRaster, [@row];
			$currow++;
			
		}
		

	}
	
	
	close $infile;
	#Raster has to be transposed as reading the ASCII leads to y-x coordinates, but MasRaster uses x-y coordinates
	# further y counting starts in lower left corner, not upper left
	my @TransposedPreRaster;
	for my $i (0..$#PreRaster){
		for my $j (0..$#{$PreRaster[$i]}){
			$TransposedPreRaster[$j][$#PreRaster - $i] = $PreRaster[$i][$j];
		}
	}

	$attributes{Raster} = MasRaster->create_from_preinitialized(\@TransposedPreRaster);
	my $self = $class->SUPER::create(%attributes);
	return $self;
}

sub reset_nodata_to {
	my $self = shift,
	my $newnodata = shift;
	my %reclasstab;
	my $oldnodatakey = sprintf("%d" ,$self->get_attribute("NoData")) ;
	#print $oldnodatakey;
	#print "\n";
	$reclasstab{$oldnodatakey} = $newnodata;
	#print $oldnodatakey. "\t $reclasstab{$oldnodatakey} \n";
	$self->{Raster}->reclass(\%reclasstab);
	$self->set_attribute("NoData",$newnodata);
}
sub round_all {
	my $self = shift;
	my $newnodata = sprintf("%d" ,$self->get_attribute("NoData")) ;
	$self->{Raster}->round();
	$self->set_attribute("NoData",$newnodata);
}

sub reclassify {
	my $self = shift;
	my $reclasstab = shift;
	my $others_to = shift;
	if (defined $others_to and $others_to eq "NoData") {
		my $missing = $self->get_attribute("NoData");
		$self->{Raster}->reclass($reclasstab,$missing);	
	}
	elsif (defined $others_to) {
		$self->{Raster}->reclass($reclasstab,$others_to);	
	}
	else {
		$self->{Raster}->reclass($reclasstab);	
	}
	$self->refresh_polygonraster() if (exists $self->{PolygonRaster});
	$self->update_valuetable() if (exists $self->{ValueTable});
}

sub random_reclassify {
	my $self = shift;
	my $reclasstab = shift;
	my $others_to = shift;
	if (defined $others_to and $others_to eq "NoData") {
		my $missing = $self->get_attribute("NoData");
		$self->{Raster}->random_reclass($reclasstab,$missing);	
	}
	elsif (defined $others_to) {
		$self->{Raster}->random_reclass($reclasstab,$others_to);	
	}
	else {
		$self->{Raster}->random_reclass($reclasstab);	
	}
	$self->refresh_polygonraster() if (exists $self->{PolygonRaster});
	$self->update_valuetable() if (exists $self->{ValueTable});
}
sub create_copy_and_reclassify {
      my $self = shift;
      my $new = $self->create_copy(); 
      $new->reclassify(@_);
      return $new;
}

sub overlay {
	my $self = shift;
	croak "EnhancedGrid->overlay cannot be called as a class method!" unless (ref $self);
	my $operation = shift;
	my $ogrid = shift;
	my $treat_missing_as_normal = shift;	
	#checks
	
	# check if both grids have the same extent
	for my $i (qw/ncols nrows xllcorner yllcorner cellsize/) {
		my $r = round_to ($ogrid->get_attribute($i), 0.01);
		my $s = round_to ( $self->get_attribute($i), 0.01);
		croak "Error when overlaying: grids do not have the same extent!\n $i is different!\n Overlay: $r \t Self: $s \n" unless ($s == $r);
	}
	# read missing values in order to 
	my $missing1 = undef;
	my $missing2 = undef;
	unless ($treat_missing_as_normal) {
		$missing1 = $self->get_attribute("NoData");
		$missing2 = $ogrid->get_attribute("NoData");
	}
	
	# call the MasRaster overlay method
	my $ogridref = $ogrid->get_attribute("Raster");
	$self->{Raster}->overlay($operation, $ogridref, $missing1, $missing2);

}
sub create_copy_and_overlay {
	my $self = shift;
	my $new = $self->create_copy(); 
	$new->overlay(@_);
	return $new;
}

sub report_cells_equal_to {
	my $self = shift;
	$self->{Raster}->report_cells_equal_to(@_);
}

sub draw_rectangle {
	my $self = shift;
	$self->{Raster}->draw_rectangle(@_);
}
sub out_of_x_bounds {
	my $self = shift;
	my @values = @_;
	$self->{Raster}->out_of_x_bounds(@values);
}
sub out_of_y_bounds {
	my $self = shift;
	my @values = @_;
	$self->{Raster}->out_of_y_bounds(@values);
}

sub tabulate_values_by_grid {
	my $self = shift;
	my $extractmap = shift;
	my $outhash = {};
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my $em_ncols	= $extractmap->get_attribute("ncols");
	my $em_nrows	= $extractmap->get_attribute("nrows");	
	print "Error in tabulate values by grid: Category map does not fit Value map\n" and return undef unless ($ncols == $em_ncols and $nrows == $em_nrows);

	for my $x (0..$ncols -1) {
		for my $y (0..$nrows -1) {
			my $value = $self->get_cell($x,$y);
			my $category = $extractmap->get_cell($x,$y);
			$outhash->{$category}->{$value} = 0 unless (exists $outhash->{$category}->{$value});
			$outhash->{$category}->{$value} += 1;
		 }
	}
	return $outhash;
}
sub sum_values_by_grid {
	my $self = shift;
	my $extractmap = shift;
	my $outhash = {};
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my $em_ncols	= $extractmap->get_attribute("ncols");
	my $em_nrows	= $extractmap->get_attribute("nrows");	
	print "Error in tabulate values by grid: Category map does not fit Value map\n" and return undef unless ($ncols == $em_ncols and $nrows == $em_nrows);

	for my $x (0..$ncols -1) {
		for my $y (0..$nrows -1) {
			my $value = $self->get_cell($x,$y);
			my $category = $extractmap->get_cell($x,$y);
			$outhash->{$category} += $value;
			
		 }
	}
	return $outhash;
}
sub set_missing_if {
      my ($self,$condition) = @_;
      $self->{Raster}->set_missing_if($condition, $self->{NoData});
}

sub list_cells_with_value {
	my ($self, $value) = @_ ;
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my @list;
	for my $x (0..$ncols -1) {
		for my $y (0..$nrows -1) {
			 push @list, {( x => $x, "y" => $y, value => $value) } if ($self->get_cell($x,$y) == $value);
		 }
	}
	return \@list;
}
sub list_cells_by_value {
	my ($self) = @_ ;
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my %list;
	for my $x (0..$ncols -1) {
		for my $y (0..$nrows -1) {	
			 my $value = $self->get_cell($x,$y);
			 push @{$list{$value}}, {( x => $x, "y" => $y)} if ($value != $self->get_attribute("NoData") );
		 }
	}
	return \%list;
}
sub list_cells {
	my ($self) = @_ ;
	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my @list;
	for my $x (0..$ncols -1) {
		for my $y (0..$nrows -1) {
			 push @list, {( x => $x, "y" => $y, value => $self->get_cell($x,$y)) } ;
		 }
	}
	return \@list;
}

sub get_header_as_string {
	my $self = shift;

	my $ncols	= $self->get_attribute("ncols");
	my $nrows	= $self->get_attribute("nrows");
	my $xllcorner	= $self->get_attribute("xllcorner");
	my $yllcorner	= $self->get_attribute("yllcorner");
	my $cellsize	= $self->get_attribute("cellsize");
	my $NODATA_value	= $self->get_attribute("NoData");
	
	return <<HD2;
ncols		$ncols
nrows		$nrows
xllcorner	$xllcorner
yllcorner	$yllcorner
cellsize	$cellsize
NODATA_value	$NODATA_value
HD2
}
1;
