[R-sig-Geo] R to Python - rpy - definition of window pattern, related with 'spatstat library - help needed

[Roger Bivand]

On Mon, 4 Feb 2008, volkan kepoglu wrote:

> I am new in R and Python. I write a couple of lines. The code does the
> followings; Read point shp file, Compute a kernel smoothed intensity
> function from a point pattern; density function, and convert pixel
> image to spatialgriddataframe and export spatialgriddataframe to tif
> raster format.
>
> the code is running in R. I am trying to use the code in python, but
> could not. Could you help me?

Cross-posting isn't a brilliant idea, by the way. I'm only replying to the 
R-sig-geo list; please summarise separately to other lists that you 
included to close the threads in a responsible way.

Do as much as possible on the R side, do not try to cross the R/python 
interface with anything other than very simple objects (scalar, numeric 
vector), unless you really enjoy debugging (and since you are writing to 
multiple lists, you probably don't enjoy debugging). Just keep things very 
simple!

>
> my original R code;

was not optimal anyway. Please use a canned file for your example:

library(maptools)
shpPoint <- readShapePoints(system.file("shapes/baltim.shp",
   package="maptools")[1])
bbox(shpPoint)
library(spatstat)
MYpppFormat <- as(shpPoint["STATION"], "ppp")
summary(MYpppFormat)
# note that window matches bounding box
z <- density.ppp(MYpppFormat)
# consider setting the bandwidth
z_sgdf <- as(z, "SpatialGridDataFrame")
outfilename <- paste(tempdir(), "exmPnt_density.tif",
   sep=.Platform$file.sep)
library(rgdal)
writeGDAL(z_sgdf, outfilename, drivername = "GTiff")

The only things that need to go over the interface are the input shapefile 
and output GTiff names, and possibly the bandwidth sigma. Try first as an 
os.system() call to a script running R, and move towards integrating 
tighter things that can be moved across the interface without falling 
over. I have a feeling that

shpPoint = r.readShapePoints ("exmPnt.shp")

has lost its class when it goes back across. I'm afraid you may get into a 
lot of r.assign() and r.get(). Best write a wrapper function taking the 
things you can pass, and leave everything else on the R side.

If anyone has implemented reflected sp objects in python, please say so!

Roger

>
> library(maptools)
> shpPoint <- readShapePoints ("exmPnt.shp")
> library(spatstat)
> w <- as.owin(c((shpPoint at bbox["coords.x1","min"]),
> (shpPoint at bbox["coords.x1","max"]),
> (shpPoint at bbox["coords.x2","min"]),
> (shpPoint at bbox["coords.x2","max"])))
> MYpppFormat <- as.ppp(shpPoint at coords, w)
> z <- density.ppp(MYpppFormat, edge=TRUE)
> plot(z, main="Kernel smoothed intensity of point pattern")
> plot(shpPoint, add=TRUE)
> z_sgdf <- as.SpatialGridDataFrame.im(z)
> outfilename <- tempfile(pattern="file", tmpdir = tempdir())
> writeGDAL(z_sgdf, outfilename, drivername = "GTiff")
> file.rename (outfilename, "exmPnt_density.tif")
>
> conversion of R code to Python with Rpy;
>
> from rpy import *
> r.library('maptools')
> shpPoint = r.readShapePoints ("exmPnt.shp")
> r.library('spatstat')
>
> # rest of them could not be converted to Rpy
> w = r.as_owin(r.c((r.bbox(shpPoint)["coords.x1","min"]),
> (r.bbox(shpPoint)["coords.x1","max"]),
> (r.bbox(shpPoint)["coords.x2","min"]),
> (r.bbox(shpPoint)["coords.x2","max"])))
> MYpppFormat = r.as_ppp(shpPoint at coords, w)
>
> # possible conversion, not tried.
> z = r.density_ppp(MYpppFormat, edge=TRUE)
> r.plot(z, main="Kernel smoothed intensity of point pattern")
> r.plot(shpPoint, add=TRUE)
> z_sgdf = r.as_SpatialGridDataFrame_im(z)
> r.library('rgdal')
> outfilename = r.tempfile(pattern="file", tmpdir = r.tempdir())
> r.writeGDAL(z_sgdf, outfilename, drivername = "GTiff")
> r.file.rename (outfilename, "exmPnt_density.tif")
>
> IN R, working;
>> w <- as.owin(c((shpPoint at bbox["coords.x1","min"]), (shpPoint at bbox["coords.x1","max"]), (shpPoint at bbox["coords.x2","min"]), (shpPoint at bbox["coords.x2","max"])))
>> w
> window: rectangle = [575490, 619430] x [4890400, 4934800] units
>
> IN PYTHON, giving the following error;
>>>> set_default_mode(NO_CONVERSION)
>>>> w = r.as_owin(r.cbind((r.bbox(shpPoint)["coords.x1","min"]),
> (r.bbox(shpPoint)["coords.x1","max"]),
> (r.bbox(shpPoint)["coords.x2","min"]),
> (r.bbox(shpPoint)["coords.x2","max"])))
> Traceback (most recent call last):
>  File "<stdin>", line 1, in <module>
> TypeError: sequence index must be integer, not 'tuple'
>>>> set_default_mode(BASIC_CONVERSION)
>>>> w = r.as_owin(r.cbind((r.bbox(shpPoint)["coords.x1","min"]),
> (r.bbox(shpPoint)["coords.x1","max"]),
> (r.bbox(shpPoint)["coords.x2","min"]),
> (r.bbox(shpPoint)["coords.x2","max"])))
> Traceback (most recent call last):
>  File "<stdin>", line 1, in <module>
> IndexError: each subindex must be either a slice, an integer,
> Ellipsis, or NewAxis
>>>>
>
> What should I do?
>
> Thanks,
> Volkan Kepoglu
> PHD Candidate,
> Department of GGIT in METU,
> Ankara, Turkey.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no