[STATSGRASS] Scripting R from outside of GRASS

[Thomas Adams]

Bill,

Attached is a bash shell script that:

(1) sets some environment variables
(2) reads some command line arguments
(3) does a GRASS r.mapcalc calculation
(4) and loops through a list of county FIPS codes (for GRASS raster masks)
(5) runs a batch R script (which just gets some summary stats)

Hope this helps.

Regards,
Tom



William Hudspeth wrote:

>Hello,
>
>I am building an online GIS application using Minnesota MapServer that
>references GRASS rasters in the map file. I have successfully displayed
>GRASS rasters, and used python to script processing of GRASS raster
>layers. I would also like to script the statistical analysis and display
>of graphs using R, but am having trouble figuring out how to run R on my
>GRASS location from outside of GRASS. I can't run R in the same script
>that I use to process GRASS rasters. Can anyone help to guide me in how
>I go about doing this. My current work flow is something like this:
>
>python_script_file.py --> calls script.sh file to load a GRASS location
>and process rasters
>
>python_script_file.py  --> call another script.sh file to generate some
>summary statistics and generate PNG graphs on the same GRASS location
>
>Any help greatly appreciated.
>
>Thanks, Bill
>
>_______________________________________________
>statsgrass mailing list
>statsgrass at grass.itc.it
>http://grass.itc.it/mailman/listinfo/statsgrass
>  
>

-- 
Thomas E Adams
National Weather Service
Ohio River Forecast Center
1901 South State Route 134
Wilmington, OH 45177

EMAIL:	thomas.adams at noaa.gov

VOICE:	937-383-0528
FAX:	937-383-0033

-------------- next part --------------
#!/bin/bash

###############################################################
# *************       ffg.xmrg.compare            *************


###############################################################
# Set GRASS Environment Variables
###############################################################
GISDBASE=/grass/data
MAPSET=oper
LOCATION_NAME=ohrfc-precip

echo "LOCATION_NAME: $LOCATION_NAME"	 > $HOME/.grassrc5
echo "MAPSET: $MAPSET"			>> $HOME/.grassrc5
echo "DIGITIZER: none"			>> $HOME/.grassrc5
echo "GISDBASE: $GISDBASE"		>> $HOME/.grassrc5
echo "GRASS_GUI: text"			>> $HOME/.grassrc5

export GISBASE=/grass/grass5sources/grass53_exp_2004_04_24/dist.i686-pc-linux-gnu
export GISRC=$HOME/.grassrc5
export PATH=$PATH:$GISBASE/bin:$PATH:$GISBASE/scripts


###############################################################
# Script to calculate the difference between 1-hr FFG & 1-hr
# XMRG grids by County (FIPS code) using GRASS & R:
# 
# 	diff = XMRG - factor*FFG
# 
#  where, factor scales the FFG to some threshold; default is factor = 1.0
# 
#  Output consists of an ASCII text file summarizing the gridded
#  differences by county FIPS code.
###############################################################
# The 'counties' variable is the full path name of the file 
# containing the list of county FIPS codes
# 
# Script converts the xmrg & ffg files to an Lambert Conformal Conic 
# Projection grid. An intermediate projection to a Lat/Long
# projection is necessary. The data is read into GRASS initially
# as a sites file
###############################################################
	
	today=$(date +%Y%m%d)
	
	counties=$1
	xmrg=$2
	ffg=$3
	factor=$4
	
	echo
	echo "=========  Calculate precipitation excess grid =  $xmrg - $factor * $ffg  ==========="
	echo
	
	cd $GISDBASE
	# Initialize the summary statistics files
	#cp /dev/null ffg.summary.$xmrg
	cp /dev/null ffg.summary
	
	# Remove any mask that may be set
	g.remove rast=MASK
	
	# Import the CPC data into GRASS
	r.mapcalc precip.excess.$xmrg=$xmrg-$factor*$ffg
	

	echo
	echo "=========  Begin processing for each Forecast Group ==========="
	echo


for county in $(<$counties); do

	echo "==========  Set mask for County $county =========="
	g.remove rast=MASK
	g.copy rast=$county,MASK
	
	export R_COUNTY_ID=$county
	export R_YEAR=$today
	export R_INMAP=precip.excess.$xmrg
	R CMD BATCH R.ffg.batch

	echo
	echo "=========  Completed for county $county  ========="
	echo

	# Cleanup -- Remove temporary GRASS files
	$GISBASE/etc/clean_temp

done

mv ffg.summary $xmrg.ffg.summary

echo
echo "==========  All counties completed  =========="
echo