If you have not already done so, you’ll need to start R and load the
spatstat
package by typing
library(spatstat)
We will study a dataset that records the locations of Ponderosa Pine
trees (Pinus ponderosa) in a study region in the Klamath
National Forest in northern California. The data are included with
spatstat
as the dataset ponderosa
.
assign the data to a shorter name, like X
or
P
;
To assign the data to X
we simply write:
X <- ponderosa
plot the data;
To plot the data we do the following:
plot(X)
By the S3 method dispatch, this calls the plot.ppp
function. The chars
argument indicate that the point type
should be be periods (“.
”).
find out how many trees are recorded;
Both npoints
and print.ppp
displays the
number of recorded trees:
npoints(X)
## [1] 108
print(X)
## Planar point pattern: 108 points
## window: rectangle = [0, 120] x [0, 120] metres
I.e. there are 108 trees in the dataset.
find the dimensions of the study region;
The dimensions of the observation window can be seen above. Alternatively, it can be directly assesed via
Window(X)
## window: rectangle = [0, 120] x [0, 120] metres
obtain an estimate of the average intensity of trees (number of trees per unit area).
The average intensity can be computed via
intensity.ppp
intensity(X)
## [1] 0.0075
or the more expansive summary.ppp
:
summary(X)
## Planar point pattern: 108 points
## Average intensity 0.0075 points per square metre
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 metres
##
## Window: rectangle = [0, 120] x [0, 120] metres
## Window area = 14400 square metres
## Unit of length: 1 metre
The Ponderosa data, continued:
When you type plot(ponderosa)
, the command that is
actually executed is plot.ppp
, the plot method for point
patterns. Read the help file for the function plot.ppp
, and
find out which argument to the function can be used to control the main
title for the plot;
From the documentation, the argument that controls the title is
main
as is also the case in the regular generic
plot
.
plot the Ponderosa data with the title Ponderosa Pine Trees above it;
plot(ponderosa, main = "Ponderosa Pine Trees")
from your reading of the help file, predict what will happen if we type
plot(ponderosa, chars="X", cols="green")
then check that your guess was correct;
Each point will be plotted as an green “X” and indeed:
plot(ponderosa, chars="X", cols="green")
try different values of the argument chars
, for
example, one of the integers 0 to 25, or a letter of the alphabet. (Note
the difference between chars=3
and chars="+"
,
and the difference between chars=4
and
chars="X"
).
There are subtle differences in the actual character/point types plotted. When given a string literal, the actual character is plotted as the point type.
plot(ponderosa, chars=3, cols="green")
plot(ponderosa, chars="+", cols="green")
plot(ponderosa, chars=4, cols="green")
plot(ponderosa, chars="X", cols="green")
The dataset japanesepines
contains the locations of
Japanese Black Pine trees in a study region.
Plot the japanesepines
data.
We use the generic plot
function which is dispatched to
plot.ppp
:
plot(japanesepines)
What is the average intensity (the average number of points per unit area?
The average intensity can be computed via
intensity.ppp
:
intensity(japanesepines)
## [1] 65
Using density.ppp
, compute a kernel estimate of the
spatially-varying intensity function for the Japanese pines data, using
a Gaussian kernel with standard deviation \(\sigma=0.1\) units, and store the estimated
intensity in an object D
say.
From the documentation (?density.ppp
) we see that the
following will work:
D <- density(japanesepines, sigma = 0.1)
Plot a colour image of the kernel estimate D
.
The plotting of the colour image is automatically done by dispatched
call to the plot.im
method by calling plot
on
the im
object.
plot(D, main = "")
Most plotting commands will accept the argument
add=TRUE
and interpret it to mean that the plot should be
drawn over the existing display, without clearing the screen beforehand.
Use this to plot a colour image of the kernel estimate D
with the original Japanese Pines data superimposed.
We can use the add = TRUE
functionality of the plotting
methods.
plot(D, main = "")
plot(japanesepines, add = TRUE, cols = "white", cex = 0.5, pch = 16)