[GRASS5] Re: [GRASSLIST:5971] r.area ?

[Jachym Cepicky]

what about this (attachement)

Jáchym

On Thu, Mar 03, 2005 at 02:39:05PM +1300, Hamish wrote:
> > is there some raster module, which would assign new values to raster,
> > based on area size of grouped cells? Something like the Area module in
> > Idrisi Kilimanjaro...
> 
> maybe with the r.to.vect module then vector command magic? (starting
> with d.what.vect)
> 
> or use r.mapcalc if(value < this && value > that,1,0)
> then 'r.stats -a' to get an area total.
> 
> maybe r.reclass + 'r.stats -a'?
> 
> 
> 
> Hamish

-- 
Jachym Cepicky
e-mail: jachym.cepicky at centrum.cz
URL: http://les-ejk.cz
GPG: http://www.fle.czu.cz/~jachym/gnupg_public_key/
-------------- next part --------------
#!/usr/bin/perl -w
use strict;

# #    LICENCE
# #    This script creates new map, which category values will represent 
# #    area of grouped cells
# #    Copyright (C) 2005 Jachym Cepicky
# #
# #    This program is free software; you can redistribute it and/or modify
# #    it under the terms of the GNU General Public License as published by
# #    the Free Software Foundation; either version 2 of the License, or
# #    (at your option) any later version.
# #
# #    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 the
# #    GNU General Public License for more details.
# #
# #    You should 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

# autor Jachym Cepicky  jachym.cepicky [at] centrum [dot] cz
# URL: http://les-ejk.cz
# 


unless($ENV{GISBASE}){
    print STDERR "\n\nYou have to run GRASS first!\n\n";
    exit;
}


my @arg = &load_arguments();

if ($arg[0] eq '0') {
    &print_help();
}

my @mapcalc_str = &reclass();

#print $mapcalc_str;
open (TMP,">/tmp/grass.r.area.txt") or die "Could not open file /tmp/grass.r.area.txt: $!\n";

print TMP @mapcalc_str;
close TMP;

system("cat /tmp/grass.r.area.txt|r.reclass in=$arg[0] out=$arg[1]; rm /tmp/grass.r.area.txt");


sub reclass 
{
    my $line = ""; # line with category from r.report
    my $cat = 0; # category number
    my $area = 0; # category area
    my @mapcalc_str = ();
    
    foreach $line (`r.report -heqnN map=$arg[0] units=$arg[2]`) {
        next unless ($line =~ m/(\. )+/);
        
#|   13| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .|      1.00|
        # searching for category
        $cat = $line;
        $cat =~ s/^\|\s*([0-9]+)\|.*$/$1/;
        chomp $cat;
        # searching for area
        $line =~ m/\. ?\|\s*(([0-9]+,)?[0-9]+(\.[0-9]+))\|$/;
        $area = $1;
        $area =~ s/,//g;

        #print "$area\n";
        push @mapcalc_str, "$cat = $area\n";
    }
    push @mapcalc_str, "end\n";
    return @mapcalc_str;
}


sub load_arguments{

    my $in="";
    my $out="";
    my $units="";
    
    # arguments of command line
    foreach(@ARGV){
        if(m/help/i){
         
            return 0;
        
        } elsif(m/in/i){ s/^.*=//;$in = $_;
        } elsif(m/out/i){ s/^.*=//;$out = $_;
        } elsif(m/un/i){ s/^.*=//;$units = $_;
        }

     }

     # testing of parameters
    if ($in eq ""){ 
        print STDERR "\nERROR: Parameter input not defined\n\n"; 
        return 0;
    }
    if ($out eq ""){ 
        print STDERR "\nERROR: Parameter output not defined\n\n"; 
        return 0;
    }
    if ($units eq ""){ 
        print STDERR "\nERROR: Parameter units not defined\n\n"; 
        return 0;
    }
    # all arguments back
    
    return ($in, $out, $units);
}



sub print_help 
{
       print qq(
\nDescription:\n
This skript creates new map, in which groups will have area size as cat. label


Usage:
 r.area in=name out=name units=name

Parameters:
 in      Name of input vectorfile, containing points
 output  Name of resulting vector file
 units   mi(les),me(ters),k(ilometers),a(cres),h(ectares),c(ell_counts),p(ercent_cover)

);

exit;
}