[STATSGRASS] grass display using R to stretch

[Agustin Lobo]

HI!,

Not really stats, but as this merges grass and R and
I'm not in the other grass lists, I send it here. Please
forward to the other lists if you feel it's worthwhile.

I've made 2 simple cshell scripts to display
grey and color images with d.rast and d.rgb
but using R to apply a 2% stretching to the
histogram. In this version the scripts do not
keep the new colr files, just use them for
display and then write back the original
ones. But modifying this would be trivial.

Also, perhaps 2% at each tail of the
histogram is too much, 1 or 1.5% might be
better.

grstretch.scr raster.in
- calculates the stretched LUT using R,
- displays (using xv) an R graph with histogram info and
the 2% and 98% values of the histogram
used to calculate the stretched LUT
- displays the stretched image using d.rast
- reverts the colr table to its original file

gr.str.rgb.scr rasterR rasterG rasterB
- calculates the stretched LUT using R for each
raster.
- displays the stretched image using d.rgb
- reverts the colr tables to their original files

Hope they are useful and other people make 
them better:

GREY SCALE IMAGERY:

grstretch.scr:

#!/bin/csh -f
#set echo on
#set verbose on

if ($#argv == 0) then
        echo "### USO: grstretch.scr raster.in ###"
        exit
endif

r.support -r $argv[1]
setenv infile /cell_misc/$argv[1]/histogram

#echo "infile <- system('echo $infile",intern=T)' > provi.in

echo "gr.stretch()" > provi.in
echo "q()" >> provi.in

R -q --save < provi.in

cp $LOCATION/colr/$argv[1] $LOCATION/colr/$argv[1].backup
cp outfile $LOCATION/colr/$argv[1]
d.rast $argv[1]
rm $LOCATION/colr/$argv[1]
mv $LOCATION/colr/$argv[1].backup $LOCATION/colr/$argv[1]
rm outfile
xv -rotate -90 Rplots.ps #si pongo &, borra antes de poder ejecutar xv
rm Rplots.ps
rm provi.in
unsetenv infile

COLOR IMAGERY

gr.str.rgb.scr:

#!/bin/csh -f
#set echo on
#set verbose on

if ($#argv == 0) then
        echo "### USO: gr.str.rgb.scr rasterR rasterG rasterB ###"
        exit
endif

##RED##
r.support -r $argv[1]
setenv infile /cell_misc/$argv[1]/histogram
echo "res <- gr.stretch()" > provi.in
echo "readline()" >> provi.in
R --no-save -q < provi.in
cp $LOCATION/colr/$argv[1] $LOCATION/colr/$argv[1].backup
cp outfile $LOCATION/colr/$argv[1]
rm outfile
unsetenv infile
##GREEN##
r.support -r $argv[2]
setenv infile /cell_misc/$argv[2]/histogram
echo "res <- gr.stretch()" > provi.in
echo "readline()" >> provi.in
R --no-save -q < provi.in
cp $LOCATION/colr/$argv[2] $LOCATION/colr/$argv[2].backup
cp outfile $LOCATION/colr/$argv[2]
rm outfile
unsetenv infile
##BLUE##
r.support -r $argv[3]
setenv infile /cell_misc/$argv[3]/histogram
echo "res <- gr.stretch()" > provi.in
R --no-save -q < provi.in
cp $LOCATION/colr/$argv[3] $LOCATION/colr/$argv[3].backup
cp outfile $LOCATION/colr/$argv[3]
rm outfile
unsetenv infile

d.rgb red=$argv[1] green=$argv[2] blue=$argv[3]

rm $LOCATION/colr/$argv[1]
mv $LOCATION/colr/$argv[1].backup $LOCATION/colr/$argv[1]
rm $LOCATION/colr/$argv[2]
mv $LOCATION/colr/$argv[2].backup $LOCATION/colr/$argv[2]
rm $LOCATION/colr/$argv[3]
mv $LOCATION/colr/$argv[3].backup $LOCATION/colr/$argv[3]
rm provi.in


THE R FUNCTION IS JUST:

"gr.stretch" <- function()
{
        layout(matrix(c(1,1,2,3), ncol=2, byrow = TRUE))
        infile.nom <-system("echo $infile",intern=T)
        print(infile.nom)
        infile <- paste(system("echo
$LOCATION",intern=T),infile.nom,sep="")
        histog <- round(matrix(scan(infile,sep=":"),byrow=T,ncol=2))
        cumfreq <- cumsum(histog[,2]*100/sum(histog[,2]))
        plot(histog[,1],histog[,2]*100/sum(histog[,2]))
        cumfreq[cumfreq<2.0] <- 200
        minim.val <- histog[cumfreq==min(cumfreq),2]*100/sum(histog[,2])
        minim <- histog[cumfreq==min(cumfreq),1]
        cumfreq <- cumsum(histog[,2]*100/sum(histog[,2]))
        cumfreq[cumfreq>98.0] <- -100
        maxim.val <- histog[cumfreq==max(cumfreq),2]*100/sum(histog[,2])
        maxim <- histog[cumfreq==max(cumfreq),1]
        title(paste("2% MIN=",minim,"98% MAX=",maxim))
        points (c(minim,maxim),c(minim.val,maxim.val),col="red")
        cumfreq <- cumsum(histog[,2]*100/sum(histog[,2]))
        plot(histog[,1],cumfreq)
        minim.val <- cumfreq[histog[,1]==minim]
        maxim.val <- cumfreq[histog[,1]==maxim]
        points (c(minim,maxim),c(minim.val,maxim.val),col="red")
        colortable <- cbind(0:255,rep(0,256))
        colortable[colortable[,1]<minim,2] <- 0
        colortable[colortable[,1]>maxim,2] <- 255
        colortable[(minim+1):(maxim+1),2] <-
rescale(minim:maxim,minim,maxim,0,255)
        plot(colortable[,1],colortable[,2])
        points(colortable[c((minim+1),(maxim+1)),1],
colortable[c((minim+1),(maxim+1)),2],col="red")
        write("% 0 255",file="outfile")

write(paste(colortable[,1],round(colortable[,2]),sep=":"),file="outfile",ncol=1,append=T)
}



Agus

Dr. Agustin Lobo
Instituto de Ciencias de la Tierra (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona SPAIN
tel 34 93409 5410
fax 34 93411 0012
alobo at ija.csic.es