Replicating the behaviour of the plot function applied to spatstat objects

[APoinas]

A lot of objects in spatstat have their own behaviour when you run "plot(name_of_the_object)" but I'm having trouble replicating this behaviour in python. For example, here is a code in R doing a rank envelope test to test if some data is generated from a Poisson PP.

library(spatstat)
rho0 = 100
alpha0 = 0.05
S_length = 2
# Generating data from a very repulsive DPP.
S = simulate(dppBessel(lambda=rho0, alpha=alpha0, sigma=0, d=2), 1, W=owin(c(0,S_length), c(0,S_length)))
# Doing an envelope test using Ripley's K function.
Env_test = envelope(S, Kest, nsim=199)
# Plotting the result of the test.
plot(Env_test)

The resulting plot is then

But, when I tried to do the same thing with spatstat-interface using the following code:

import numpy as np
import pandas as pd

import rpy2.robjects as robjects
from rpy2.robjects import numpy2ri
numpy2ri.activate()

from spatstat_interface.interface import SpatstatInterface
from spatstat_interface.utils import convert_r_df_to_pandas_df

spatstat = SpatstatInterface(update=True)
spatstat.import_package("core", "geom", update=False)

rho = 100
alpha = 0.05
sigma = 0
d = 2
params = {"lambda": rho, "alpha": alpha, "sigma": sigma, "d": d}
my_dpp = spatstat.core.dppBessel(**params)
bound = robjects.FloatVector([0, 2])
window = spatstat.geom.owin(xrange=bound, yrange=bound)
# Generating data from a very repulsive DPP.
rsample = spatstat.core.simulate_dppm(my_dpp, W=window)
# Doing an envelope test using Ripley's K function.
Envelope = spatstat.core.envelope(rsample, spatstat.core.Kest, nsim=199)
# Plotting the result of the test.
base = robjects.packages.importr("base")
base.plot(Envelope)

The resulting plot is then which is definitely not what we want.

Could there be a way we can replicate with spatsat-interface the behaviour of the plot function applied to spatstat objects?

[guilgautier]

Hi @APoinas! I have just activated the "Watch:All activity" and should now be warned when issues are raised 😅

Regarding your problem, have you tried the spatstat.core.plot.envelope function ?

[APoinas]

No problem! Unfortunately, when I try to call spatstat.core.plot.envelope I get AttributeError: module 'spatstat.core' has no attribute 'plot'

Moreover, the issue doesn't appear just for "envelope" type objects. For example,

library(spatstat)

Kestim = Kest(rpoispp(100)) 
plot(Kestim)

in R gives while

import numpy as np
import pandas as pd

import rpy2.robjects as robjects
from rpy2.robjects import numpy2ri
numpy2ri.activate()

from spatstat_interface.interface import SpatstatInterface
from spatstat_interface.utils import convert_r_df_to_pandas_df

spatstat = SpatstatInterface(update=True)
spatstat.import_package("core", "geom", update=False)

Kestim = spatstat.core.Kest(spatstat.core.rpoispp(100))
base = robjects.packages.importr("base")
base.plot(Kestim)

in python gives Here, the plot function is applied to an object of class "fv".

As another example consider this R code simulating a DPP and plotting it:

library(spatstat)

rho0 = 100
alpha0 = 0.05
S_length = 2
S = simulate(dppBessel(lambda=rho0, alpha=alpha0, sigma=0, d=2), 1, W=owin(c(0,S_length), c(0,S_length)))
plot(S)

gives When trying to do the same in python:

import numpy as np
import pandas as pd

import rpy2.robjects as robjects
from rpy2.robjects import numpy2ri
numpy2ri.activate()

from spatstat_interface.interface import SpatstatInterface
from spatstat_interface.utils import convert_r_df_to_pandas_df

spatstat = SpatstatInterface(update=True)
spatstat.import_package("core", "geom", update=False)

rho = 100
alpha = 0.05
sigma = 0
d = 2
params = {"lambda": rho, "alpha": alpha, "sigma": sigma, "d": d}
my_dpp = spatstat.core.dppBessel(**params)
bound = robjects.FloatVector([0, 2])
window = spatstat.geom.owin(xrange=bound, yrange=bound)
rsample = spatstat.core.simulate_dppm(my_dpp, W=window)
base.plot(rsample)

we get which is closer to the R output but still different. And in this case, the object we plot is of type "detpointprocfamily".

[guilgautier]

First let's compare what is comparable

Initial code

R code

library(spatstat)
print(packageVersion("spatstat.core"))  # 2.3.0

rho = 100
alpha = 0.05
sigma = 0
d = 2

my_dpp = spatstat.core::dppBessel(lambda=rho, alpha=alpha, sigma=sigma, d=d)
bound = c(0, 2)
window = spatstat.geom::owin(xrange=bound, yrange=bound)
rsample = spatstat.core::simulate.dppm(my_dpp, W=window)
env = spatstat.core::envelope(rsample, spatstat.core::Kest, nsim=199)

spatstat.core::plot.envelope(env)

Python code

import rpy2.robjects as robjects

from spatstat_interface.utils import convert_r_df_to_pandas_df
from spatstat_interface.interface import SpatstatInterface

spatstat = SpatstatInterface(update=True)
spatstat.import_package("core", "geom", update=True)
spatstat.core.__version__ # 2.3-0

rho = 100
alpha = 0.05
sigma = 0
d = 2

params = {"lambda": rho, "alpha": alpha, "sigma": sigma, "d": d}
my_dpp = spatstat.core.dppBessel(**params)
bound = robjects.FloatVector([0, 2])
window = spatstat.geom.owin(xrange=bound, yrange=bound)
rsample = spatstat.core.simulate_dppm(my_dpp, W=window)
env = spatstat.core.envelope(rsample, spatstat.core.Kest, nsim=199)

spatstat.core.plot_envelope(env)   # note the _ instead of the . in plot_envelope

---

In both cases I obtain

Error in NextMethod("plot", main = main) : 
  'NextMethod' called from an anonymous function

I find this very weird, but let's assume that's the normal behavior.

---

Plotting

With R

base::plot(env)

renders a nice plot as you posted above

With Python

base = robjects.packages.importr("base")
base.plot(env)

doesn't render any plot, but

R/rpy2 DataFrame (4 x 5)
lty col key label   meaning
... ... ... ... ...

where

env 

R/rpy2 DataFrame (513 x 5)
r | obs | theo | lo | hi
-- | -- | -- | -- | --
... | ... | ... | ... | ...

in a VSCode notebook. What is your IDE ?

[guilgautier]

More generally, I think the plotting has nothing to do with the spatstat-interface package as it is a very simple wrapper based on rpy2. We should probably have a look at how rpy2 handles graphics, see, e.g., the graphics section in rpy2's documentation.

[guilgautier]

Nevertheless, one can simply use pandas.DataFrame as the Python counterpart of R DataFrame Simply run the following, you should get something that looks similar to what you're looking for

from spatstat_interface.utils import convert_r_df_to_pandas_df
df = convert_r_df_to_pandas_df(env)
df.plot(x="r")

[APoinas]

Yeah, it's pretty weird that we get different results. My IDE is Spyder so it could be IDE related.

I didn't know I could call plot.envelope in R by using plotenvelope with a . Although, when I run your first python code in https://github.com/For-a-few-DPPs-more/spatstat-interface/issues/1#issuecomment-934592409 I get a different error message than you do: RRuntimeError: Error in verifyclass(X, "fv") : argument ‘X’ is not of class ‘fv’ It's pretty weird.

[guilgautier]

Although, when I run your first python code in #1 (comment) I get a different error message than you do: RRuntimeError: Error in verifyclass(X, "fv") : argument ‘X’ is not of class ‘fv’ It's pretty weird.

The reason may be that you don't work with the latest version of the packages

# R
print(packageVersion("spatstat.core"))  # 2.3.0

# Python
spatstat.core.__version__ # 2.3-0

I didn't know I could call plot.envelope in R by using plotenvelope with a .

In R it's with a ., in Python using rpy2 it's a _

[APoinas]

You're right, it seems that I have the latest version of spatstat in R but not in Python where it's version 2.2.0 that is loaded. Using "update=true" when loading the packages in python doesn't seem to change it. Although, at least that explains the different outputs.

[APoinas]

After internal discussions, it seems that it's difficult to replicate the general behaviour of the plot function in R. Although, there seems to be some workaround.

Using spatstat.core.plot_fv(env) works for plotting envelopes; Using spatstat.core.plot_fv(Kest) works for plotting Ripley's function estimates.

Also, the version of spatstat.core absolutely needs to be 2.3-0.