User Tools

Site Tools


wiki:langintergartionptot

Languages/Software data integration: GRASS and R


The main aim of this page is to show how to integrate GRASS and R sw under BASH environment.
For this we create an exercise using ~/ost4sem/exercise/basicadvgdalogr/fagussylvatica/20*43023435.tif and ~/ost4sem/exercise/basicadvgdalogr/fagussylvatica/polyfr.shp. The object of the exercise is to calculate the Moran's I index (using R) inside the polygons of the polyfr.sh shape file.

Open the files with openev to get familiar with the geodata.

Setting WORKING directory.

 DIR=~/ost4sem/grassdbnew/

Create a location in a new grass database using a referenced dataset.


Copy referenced dataset in the working directory.

mkdir $DIR
cd $DIR
cp ~/ost4sem/exercise/basic_adv_gdalogr/fagus_sylvatica/2020_TSSP_IM-ENS-A2-SP20_43023435.tif   $DIR/2020_TSSP_IM-ENS-A2-SP20_43023435.tif

Rather than using the GRASS GUI to create a location we can run a bash script (createlocation.sh) which create a location (in our case called EUlocation) using an existing dataset (in our case 2020TSSPIM-ENS-A2-SP2043023435.tif).

The create
location.sh scripts need 3 arguments:

createlocation.sh rasterfile newlocationname [GISDBASE]

bash ~/ost4sem/exercise/basic_adv_grass/create_location.sh   2020_TSSP_IM-ENS-A2-SP20_43023435.tif  EU_location ~/ost4sem/grassdbnew 

At this point the location has been created and we can start grass in the EU_location.

Setting GRASS variables for GRASS bash job

GRASS needs to know the Location and the Mapset to start. This can be done by the GRASS GUI or by a command line ( e.g. grass -text ~/ost4sem/grassdb/europe/PCEM).
GRASS saves the starting settings in a text file called /home/user/.grassrc6. You can open this file by more /home/user/.grassrc6

We can create the .grassrc6 file based on the LOCATION and MAPSET that we want enter.

echo "LOCATION_NAME: EU_location"                >  $HOME/.grassrc6
echo "MAPSET: PERMANENT"                         >> $HOME/.grassrc6
echo "DIGITIZER: none"                           >> $HOME/.grassrc6
echo "GRASS_GUI: text"                           >> $HOME/.grassrc6
echo "GISDBASE: /home/user/ost4sem/grassdbnew"   >> $HOME/.grassrc6

Now we can set/export the paths to the GRASS binaries and libraries and we can start GRASS by using GISRC=~/.grassrc6

export GISBASE=/usr/lib/grass64
export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
export LD_LIBRARY_PATH="$GISBASE/lib"
export GISRC=~/.grassrc6

Grass create a file .gislock which prevent to use the same location under different grass section.
Therefore we can create an unique .gislock file using the $$ simbol that produce a unique process ID (PID) number.

export GIS_LOCK=$$

Test if you enter correctly in GRASS by running the g.gisenv

g.gisenv

Import all tif dataset in GRASS

At this point we are ready to import the ~/ost4sem/exercise/basicadvgdalogr/fagussylvatica/2043023435.tif files in GRASS using r.in.gdal command. <code bash> for file in ~/ost4sem/exercise/basicadvgdalogr/fagus_sylvatica/2043023435.tif ; do r.in.gdal input=$file output=$(basename $file .tif) –o done </code> To keep clean the GRASS Mapset we can remove the dataset that we used to create the location. <code bash> g.remove rast=2020TSSPIM-ENS-A2-SP2043023435.tif </code>

Import poly shape file in GRASS

Import a polygon shape file of two regions into GRASS using v.in.ogr command.

v.in.ogr dsn=/home/user/ost4sem/exercise/basic_adv_gdalogr/fagus_sylvatica/poly_fr.shp  output=poly_fr  --overwrite -o

Setting region

Set a working region using poly_fr and create small datasets (using r.mask) before importing into R.

# convert vector to raster 
v.to.rast input=poly_fr output=poly_fr  use=attr  column=id  --overwrite
# setting the region
g.region vect=poly_fr@PERMANENT res=1000
 
# create a mask using the polygon ID=1 (left polygon)
r.mapcalc "poly_fr1 = if ( poly_fr == 1 , 1 , null()  ) "
 
# apply the mask (all the following operation will be based on the region and mask settings)
r.mask input=poly_fr1 -o
 
# loop through all the rasters and 'mask' them (set all masked values to NA)
for file in `g.mlist type=rast  pattern=20*IP*35  `; do 
r.mapcalc  "'$file"_fr1"'" =  "'$file'" 
done

Calculate Moran's I in R

Create a loop that calculates Moran's I (http://en.wikipedia.org/wiki/Moran's_I) to each raster in GRASS.

# get the list of rasters with names that start with "20", followed by "IP", and ending with "_fr1"
for file_fr1 in `g.mlist type=rast  pattern=20*IP*_fr1  `; do  
export file_fr1=$file_fr1
 
 
# start R and import bash variables (the file name)
 
R --vanilla -q <<EOF
 
file_fr1 = Sys.getenv('file_fr1')
library(raster)
rasterOptions(tmpdir="/tmp")
 
# load grass raster data 
 
rastD = raster( paste ("/home/user/ost4sem/grassdbnew/EU_location/PERMANENT/cellhd/",file_fr1,sep=""))
 
# calculate Moran's I
M_text = Moran(rastD, w=9)
 
#Write the results to a simple text file
write.table (M_text[1] ,paste("moran_",file_fr1,".txt",sep="" )  ,row.names = F , col.names = F)
EOF
 
done

Look at these files (which contain Moran's I for each separate raster)

 cat $DIR/moran_*_TSSP_IP-*_43023435_fr1.txt 

Extract Moran's I average using AWK

 cat $DIR/moran_*_TSSP_IP-*_43023435_fr1.txt   | awk '{ sum=sum+$1 } END { print sum/NR}  ' 

Exercise

Using the right side polygon in the polyfr calculate the Moran's I and compare the results with the Moran's index obtained in the previous calculation. Create a unique script and try to run with bash ~/ost4sem/exercise/GRASSRintegration/grassr_integration.sh

wiki/langintergartionptot.txt · Last modified: 2017/12/05 22:53 (external edit)