[GRASSLIST:6428] Re: satellite images

[Jachym Cepicky]

what about histogram stretch?

I wrote some script, which uses R to cut out all extrem values and tan it
provides the histogram strech using r.mapcalc

maybe it helps

see http://www.fle.czu.cz/~jachym/.../histogram_strech1.png
and attachement

Jáchym

On Tue, Apr 12, 2005 at 10:11:00AM +0200, Markus Neteler wrote:
> On Mon, Apr 11, 2005 at 04:52:45PM +0200, Chiara Porcu wrote:
> > I all!
> > 
> > I've imported a satellite image en grass 60, with a resolution of 0.6 m,
> > and 16 bits. the problem is that the white/black image is too dark and the
> > color one appears just as it was taken in a very foggy day. but when I open
> > them in windows it's all ok. can it be a problem of composition? is there
> > a way to adjust them on grass? I used r.in.gdal to import the images and
> > r.composite to compose them. in a first moment we worked on images of 8 bits
> > and, even if the quality wasn't good, it appeared better than now. can it
> > be a problem of bits?
> 
> To visually improve channels, you can use r.colors (col=aspect does
> sometimes a magic improvement). I assume that you are using Quickbird.
> 
> If you have all the color and panchromatic channels you can use
> i.fusion.brovey to make an image fusion.
> 
> Markus
> 

-- 
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 --------------
#!/bin/sh

############################################################################
#
# skript:	upravads.sh
# autor:	Jachym Cepicky
# popis:	upravuje zadane druzicove snimky, histogram stretch atd..
#               do aktualni lokace
# copyright:	This program is free software under the GNU General Public
#		License (>=v2). 
# pozadavky:    - sed
#               - grep
#               - R
#############################################################################

if test "$GISBASE" = ""; then
 echo "You must be in GRASS GIS to run this program." >&2
 exit 1
fi

#%Module
#% description: Provides histogram stretch
#%End
#%option
#% key: raster
#% type: string
#% description: raster file
#% required : yes
#%End
#%option
#% key: suffix
#% type: string
#% description: suffix for the result map
#% required : yes
#%End
#%option
#% key: resmult
#% type: integer
#% description: multiplicator for the raster resolution. usefull for faster R calculations, the result raster will be untouched
#% required: no
#%End

if [ "$1" != "@ARGS_PARSED@" ] ; then
  exec g.parser "$0" "$@"
fi

#####################################################################
# global variables
#####################################################################
cat="$GIS_OPT_raster"
suff="$GIS_OPT_suffix"
resol="$GIS_OPT_resmult"

#####################################################################
# KONTROLY
#####################################################################
# kontrola sed
if [ -z "`which sed`" ] ; then
   echo "ERROR: Script needs sed."
   exit 1
fi

# kontrola grep
if [ -z "`which grep`" ] ; then
   echo "ERROR: Script needs grep."
   exit 1
fi


# kontrola R
if [ -z "`which R`" ] ; then
   echo "ERROR: Script needs R (www.r-project.org)."
   exit 1
fi


#####################################################################
# FUNKCE PROGRAMU
#####################################################################

##
## Provides histogram stretch using R statitstick package
##
## Jachym Cepicky 2005
function histogram_stretch()
{
    ds=$1
    suff=$2;
    resol=$3

    # nelze delat s celym vyrezem, resolution = 2*res
    #g.region rast=$ds; 

    res=`g.region -p |grep wres:|sed -e s/.*res\:\ *//` # aktuální rozlišení
    res=`echo $res | sed -e "s/\([0-9]*\)\..*/\1/"`; # bez desetinné tečky
    if [ "$res" == "0" ]; then
        res=1;
    fi
    res=$(( ${res} * ${resol} )); # výsledné rozlišení

    #echo "Nastavuji region na ${resol}xrozlišení rastru";
    echo "Setting region to ${resol}xraster resolution";
    
    g.region res=$res   # nastavíme
    g.region -p |grep res       # pro jistotu vytiskneme

    
    #echo "Spouštím statistický program R";
    echo "Running R";
    rout=`echo "library(GRASS); g<-gmeta(interp=TRUE);
          ds<-rast.get(g, \"$ds\",c(F), interp=T);
          attach(ds);
          stats<-boxplot.stats($ds,  do.conf = TRUE, do.out = FALSE);
          stats$stats[1];
          " | R --vanilla --no-save -q --slave|grep [1]`;
          
    #echo "Rko proběhlo, seduju výsledky";
    echo "R calculations are over, I'm seding out the results";
    rout=`echo $rout | sed -e s/.*gisrc\ .1.\ //`;  # tady jsou čísla
    min=`echo $rout  |sed -e "s/\([0-9]*\).*/\1/"`; # pouze minimum
    max=`echo $rout  |sed -e "s/.*\ \([0-9]*\)$/\1/"`; # pouze maximum

    #echo "Provedu natažení histogramu: minimum: $min->0, maximum: $max->255"
    #echo "Nastavuji region"
    echo "Setting region"
    g.region rast=$ds;

    #echo "Provádím roztažení histogramu";
    echo "Providing histogram stretch: minimum: $min->0, maximum: $max->255"
    r.mapcalc $ds.$suff="eval(out=($ds-$min)*255/($max-$min),if( out>255, 255, if (out < 0, 0, out)))"

    #echo "Nastavuji barevnou paletu na grey.eq"
    echo "Setting color palet to grey.eq"
    r.colors map=$ds.$suff color=grey.eq
    return 0;
}



#####################################################################
# TELO PROGRAMU
#####################################################################
if [ "${resol}" == "" ]; then
    resol=2;
fi

if [ "${cat}" == "" ]; then
    echo "ERROR: Input raster <raster> not set"
    exit 1;
fi
if [ "${suff}" == "" ]; then
    echo "ERROR: Output suffix <suff> not set"
    exit 1;
fi



eval `g.gisenv`
: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}

histogram_stretch $cat $suff $resol;