Jesus Loves GRASS

Technical blog about GRASS/GIS,open source geoinformatics and MAPSERVER.

Thursday, August 02, 2007

Creation of a Semi-Variogram using RSOAP and Python

This is basically what is in Py Zine about RSOAP but with some differences....

Starting in the bash shell by activating the RSOAPManager : ./RSOAPManager --debug

########################################################
STARTING R SOAP Manager on localhost:9081
########################################################
Manager Port: 9081
R Server Port range: [9082 to 65536]

Setting up R Process Object...
../lib/RProcess.py:70: RuntimeWarning: tempnam is a potential security risk to your program
self._tmpdir = os.tempnam(self.path,"RSOAP") + "/"
RHOME= /usr/local/lib/R
RVERSION= 2.5.1
RVER= 2051
RUSER= /root
Loading Rpy version 2051 .. Done.
Creating the R object 'r' .. Done
Setting up SOAP Listener

########################################################
Listening for SOAP requests at http://localhost:9081
########################################################

Waiting for request...


jumping to the python shell:

>>from SOAPpy import *

Starting a new SOAP connection:
>>RSOAPServer=SOAProxy(“http://localhost:9081”) # as indicated by RSOAPManager

To check if the connection is ok
>>RSOAPServer.echo(“something”) #the server will reply to what ever we write
'something'

Requesting a new R session and port information:
>>RSessionURL=RSOAPServer.newServer() # the reply should be http://localhost:9083 if the reply is 'None' it didnt opened the session

It is also possible to keep track of the ports using the nmap in the bash shell:
>nmap -v -p9080-9090 localhost

The making the actual connection to the R session
>>RSession=SOAPProxy(RSessionURL)

checking:
>>RSession.echo(“test”) # for killing session : RSession.quit()


3 major ways to deal with the R session: function call, eval, script. The first one (call) is used to call functions and pass parameters, eval is used to send lines of commands (ex: names(ssa)<-c(“x”,”y”,”z”)) and script to process complete scripts.

For the semivariogram creation it will be used the sgeostat pack. The first thing is to load the library

>>RSession.call(“library”,”sgeostat”) #fist it is called the functions and then the parameters to pass

Using sgeostat it is necessary the following steps to create a semivariogram
-Upload data
-Create point object
-Create pair object
-Estimate the semi-variogram

For loading the data set and creating the data object (ssa):
>>RSession.eval("ssa<-read.csv(file='/home/jesus/doutoramento/db/SSA.csv',header=FALSE,sep=';')")

The eval will pass the command R.

to see the new object:
>>RSession.call(“ls”)
'ssa'

For the other sequential commands it is better to use the script option. Normally for this method it is assumed the \n as the escape char for separation, it also allows for prompts when the echo option is activated.

script_text="names(ssa)<-c('x','y','z') \n ssa.point<-point(ssa) \n ssa.pair<-pair(ssa.point) \n ssa.estvar<-est.variogram(ssa.point,ssa.pair,'z')"

Sending the script:
>>Rsession.script(script_text)

Checking:
>>RSession.call("ls")
['resultText', 'ssa', 'ssa.estvar', 'ssa.pair', 'ssa.point']

Normally to output the result I use png() to save the result as a png file, the problem is that X11() device is deactivated (it took me sometime to discover....), to see what device is being used:

>>RSession.call(“getOption”,”device”)
'postscript'

Another possibility is to use the PDF output
>>> RSession.eval("options(device='pdf')")

Using the PDF the result will be dumped using a script that opens the PDF device, plots, closes/saves
>>script_text2="pdf(file='/home/jesus/tmp/sv.pdf') \n plot(ssa.estvar) \n dev.off()"
>> RSession.script(script_text2)

The result should be the SV PDF image with the semivariogram :)

Labels: , , ,

0 Comments:

Post a Comment

<< Home