0

I have a contour map in spatstat generated from the intensity function of a point pattern X (like "location of the trees"). Each x,y coordinates in this point pattern is marked with a corresponding third vector (like "diameter of the tree"). -->cf image (of course the vertical lines representing the tree can be omitted)

I would like to display the average of the mark (diameter) in each level of the contour with different colors. Suggestions? Thanks!

saccente
  • 21
  • 1
  • Your question is bit vague and I don't know exactly what you want. Could you please supply an example with code of what you have done so far? Possibly use a built-in dataset as an example (e.g. `spruces` which has locations and numeric marks with sizes of trees). – Ege Rubak Sep 25 '15 at 11:22
  • After reading you post again I think I understand it better: You want to extract the spatial regions defined by the contours; extract the trees in each of these regions and calculate their mean size; make a filled contour plot with colour corresponding to the mean size. If that is right I think the answer is that there is no easy way to do this. What do you need the plot for? Could you use the `weights` argument of `density.ppp` to obtain a plot (possibly in combination with other plots) that conveys the information you want? – Ege Rubak Sep 25 '15 at 11:52

1 Answers1

0

You are effectively asking for a kind of nonparametric regression.

Here is a quick-and-dirty calculation using the function rhohat and demonstrated on the longleaf dataset.

First calculate the intensity estimate: Z <- density(longleaf) yielding an image Z. Next treat Z as a covariate in calls to the rhohat command:

f <- rhohat(unmark(longleaf), Z)

and

g <- rhohat(unmark(longleaf), Z, weights=marks(longleaf)).

Now take the ratio, h <- eval.fv(g/f) and plot it, plot(h). This shows the estimated average tree diameter as a function of the forest density. To apply this function h to the original contours of Z you would first convert h to a true function by H <- as.function(h) then evaluate hZ <- eval.im(H(Z)) and finally plot(hZ).

  • Thank you Adrian, it works (not that I had any doubts since the answer was coming from you...). Looking forward for the spatstat book by the way! – saccente Oct 23 '15 at 05:32