Legend in ppp plot

I´m creating a spatial map of trees based on x-y location using the spatstat package and the ppp function. I´ve been able to do this and the size of each point is based on tree diameter. What I haven´t been able to do, and I hope to find some help here, is to create a legend that allows me to show these different circles and also include the 6 tree species within the plot. Here is the code I´m using:

df <- subset(plots, plots$spp == "DF")  # Douglas-fir

dfx <- as.numeric(as.character(df$x))
dfy <- as.numeric(as.character(df$y))
dfd <- as.numeric(as.character(df$d))

dfp <- ppp(dfx, dfy, window = owin(c(0, 100), c(0, 100)),
                  unitname=c("metres","metres"), marks = dfd)
par(mar = c(2, 2, 2, 2))                                           
plot(dfp, main = "", cex = 0.8, markscale = 0.04, 
     bg = rgb(0.1,0.9,0.3,0.5), fg = "black")

I have a similar structure for every species in the plot (wh for Western Hemlock and the code is the same...)

Thanks!

Solution

First collect all the data into a single ppp object using superimpose. For example if you have Douglas Fir dfp and Western Hemlock whp then

  X <- superimpose(DF=dfp, WH=whp)

creates a point pattern X in which each point has two mark values, the species and the diameter. For convenience, change the mark column names:

  colnames(marks(X)) <- c("diameter", "species")

Next determine the plot symbol map for the diameters (without plotting anything):

  dmap <- plot(subset(X, select=diameter), do.plot=FALSE)

Next pick a set of colours for the different species, e.g.

  spec <- levels(marks(X)$species)
  scol <- c("red", "blue")
  smap <- symbolmap(inputs=spec, col=scol)

Now plot the species in separate colours using the diameter scale map:

  plot(Window(X), main="The main title")
  for(i in seq_along(spec)) {
     sy <- update(dmap, col=scol[i])
     Xi <- subset(X, species==spec[i], select=diameter)
     plot(Xi, add=TRUE, symap=sy)
  }

Finally plot the two symbol maps dmap and smap in the places you want by using plot(dmap, add=TRUE, xlim=..., ylim=...) and so on.

In future this will be automated, but it hasn't been implemented yet.