Run R script from Python

Question

I want to run an R-script from a Python script. The R-script is required for the projection of lat lon coordinates in a different coordinate system. I have considered two options to do this. In the first option I like to parse the lat and lon coordinates to the R-script which is shown below. Then finally I would like that the R-script returns the x and y back to the python script, but I can't figure out how to do this.

project<-function(lat,lon){

library(sp)
library(rgdal)

xy <- cbind(x = lon, y = lat)
S <- SpatialPoints(xy)
proj4string(S) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
Snew <- spTransform(S, CRS("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"))
x <- coordinates(Snew)[1]
y <- coordinates(Snew)[2]

return(x,y)
}

For my second option I've considered using the R-script at the bottom with the lat and lon already in it. I try to run this from python with subprocess.Popen('Rscript project.r', shell=True).wait() But this does not seem to work. It does not write an xy.txt file. If I run this from the cmd line, however, the R-script does the job. Who can help me out with one of these two options?

library(sp)
library(rgdal)

lat <- 52.29999924
lon <- 4.76999998

xy <- cbind(x = lon, y = lat)
S <- SpatialPoints(xy)
proj4string(S) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
Snew <- spTransform(S, CRS("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"))
x <- coordinates(Snew)[1]
y <- coordinates(Snew)[2]

cat(x, file="xy.txt",sep="")
cat(y,file="xy.txt",append=TRUE)

Just out of curiosity, is there a reason why you don't just use the Python GDAL binding? (pypi.python.org/pypi/GDAL) Seems like it would be a lot easier. — rici
– rici, Commented May 8, 2013 at 23:19
If it's easier then I definitely would like to try that. However, I am not familiar with using gdal in python. Is there an example code of how to use gdal and project coordinates from one coordinate system into another? — user2357183
– user2357183, Commented May 10, 2013 at 10:46
I don't think it's much different from the R binding. If you just need to transform co-ordinates, you only need the proj4 bindings; it's as simple as this: import pyproj; p1 = pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"); p2 = pyproj.Proj("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"); pyproj.transform(p1, p2, 52.3, 4.77) The lat/lon arguments can be lists, tuples, numpy arrays or scalars. See pyproj.googlecode.com/svn/trunk/docs/index.html — rici
– rici, Commented May 10, 2013 at 17:06

Paul Hiemstra · Accepted Answer · 2013-05-09 05:03:57Z

1

It looks to me you are almost there, some remarks:

x <- coordinates(Snew)[1] selects the first item, use x <- coordinates(Snew)[,1] to get the entire vector. This assumes that you pass vectors of lat and lon, and not single values. There is no reason here why we should not allow passing vectors of lat lon.
return(x,y) is not valid, R does not support returning multiple objects. In stead, just return(cbind(x,y)).
Inside a function I would use require.

In regard to running the code from Python, I'd have a look at Rpy. This allows you to call R code from within Python. Rpy does a lot of the hard work of passing R objects back and forth for you, so this is a nice option in this case.

The approach I would take in Rpy is to put project in a .R file, source that file, and call project.

edited May 9, 2013 at 5:03

answered May 8, 2013 at 17:31

Paul Hiemstra

61.2k12 gold badges146 silver badges151 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

G. Grothendieck Over a year ago

library will also not reload a loaded package.

Community · Accepted Answer · 2017-05-23 10:25:13Z

0

For your first solution -- get R to output the coordinates to the console (using print or cat or whatever it is in R), then capture that in Python (see here: Running shell command from Python and capturing the output)

That will give you a string with teh lat/lon coords in it, you just have to Parse them out using the appropraite Python functions.

edited May 23, 2017 at 10:25

CommunityBot

11 silver badge

answered May 8, 2013 at 17:26

Colin Nichols

7144 silver badges12 bronze badges

2 Comments

Paul Hiemstra Over a year ago

There is no need to manually parse the output of the R function, Rpy can do that for you much more easy.

user2357183 Over a year ago

Hi all, I'd like to thank you all for the very good answers. I think Rpy can do the job for me, however, Rpy does not seem to work with python version 2.5.4 and R version 2.14.2?

Collectives™ on Stack Overflow

Run R script from Python

2 Answers 2

1 Comment

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related