[R-sig-Geo] Kcross question

[Adrian.Baddeley at csiro.au]

In addition to Rolf's comments:

      In spatstat version 1, all computations are performed using the Euclidean distance. If your data were originally recorded as latitude, longitude
      coordinates then ideally it would be better to measure the great circle (geodesic) distance between points, and to compute the Kcross function (etc)
      using these distances; but we have not yet released any code in spatstat to do this. We are working on it. 

      In the meantime the best strategy is to project the lat, long coordinates using a suitable projection (see the packages spproj and maptools) 
      and then to use spatstat to calculate Kcross etc.

      Adrian Baddeley
       
________________________________________
From: Rolf Turner [r.turner at auckland.ac.nz]
Sent: Sunday, 11 March 2012 12:16 PM
To: Ashly Townsen
Cc: r-sig-geo at r-project.org; Baddeley, Adrian (CMIS, Floreat)
Subject: Re: [R-sig-Geo] Kcross question

(1) At the present time "spatstat" has no facilities for for projecting one
coordinate system onto another.  There are tools, in the "sp" package,
for doing such projection, but I have no experience in using them.  Others
on this list will be able to help you.

(2) In respect of the "r" and "breaks" argument:   The help says not to mess
with "breaks" in normal circumstances. If you want to ``use your own increments''
just specify the "r" argument, making sure that the entries of "r" are non-negative,
and strictly increasing.  And the first entry of "r" is 0.

But *WHY* do you want to specify ``your own increments''?  You are advised
not to do so, and there is no indication in what you have written as how specifying
your own "r" vector would make the task of analysing your data any easier.

(3) If you want ``output that is in kilometres or metres'' then you must specify
the coordinates of your pattern in kilometres or metres.  That is, you must
project your long and lat coordinates onto a flat surface on which distance
is measured in the required units.  I *think* that ``UTM'' is the system you are,
or may be, after but I may be talking through my hat here.  Others may advise
you more wisely.

Note that if you get your coordinates set up as, e.g., kilometres, you may wish
to specify the window for your pattern using argument unitname="kilometres".
The default x-axis label for your envelope plot will then include "kilometres".

(4) It is not clear to me what you are really getting out of using the Kcross()
function anyway.  If I were analysing such data I would start by forming a
distance map for each of the patterns of battle clusters and then use that
distance map as a covariate in a model to predict the intensity of the peace
keeping operation pattern.  Using the rhohat() function could be revealing.

    cheers,

        Rolf Turner

On 11/03/12 11:58, Ashly Townsen wrote:
Hi,

I am creating a bivariate K analysis with the spatstat package.  The code I used is below, but perhaps a little explanation is preferable first.  I am trying to analyze the location of peacekeeping operations relative to battle clusters in countries experiencing a civil war (I am a poli. sci. graduate student).  The analysis should tell me if peacekeeping locations are located near different types of battle clusters (this example is low, med, high deaths) and I have created an output graph that should illustrate how each compares to the other.  The data is all GIS (long, lat) coded and I have used the country's border as the bound.  The main question I am struggling over is related to the r (radius?) value in the kcross command.  Is there a way to force the r to use my desired increments (I know the help files says this is possible)?  I would prefer an output that is in kilometers or meters.  If not, how do I determine the value of r in an unit measurement?  I have seen a few outputs with x axes labels reading r= .5 meters but I can't figure out how to do this myself.  I have tried to use the breaks argument but I can't seem to get the function work.  Is there a pre-step I need to conduct to establish my GIS data as a certain projection or tell R how my data is coded?   I have also attached the output if you wish to take a look and let me know if you need additional information.

Thanks in advance,


Ashly Adam Townsen
University of Illinois
Department of Political Science



## CREATE BOUND (FOR DRC)
drc <- readShapePoly("drc2", repair=TRUE)
bound.drc <- as.owin(drc, "owin")

## BRING BACK AS MULTITYPE PPP
data.ppp <- scanpp("DRC_BD.txt", window=bound.drc, multitype=TRUE)

#### YEAR 2004 ####

## LOW BATTLE DEATHS
bik.low.drc.sim.2004 <- envelope(data.ppp, correction="best", Kcross, i="490LowBD2004", j="490pko2004", nsim=100)

## MED BATTLE DEATHS
bik.med.drc.sim.2004 <- envelope(data.ppp, correction="best", Kcross, i="490MedBD2004", j="490pko2004", nsim=100)

## HIGH BATTLE DEATHS
bik.high.drc.sim.2004 <- envelope(data.ppp, correction="best", Kcross, i="490HighBD2004", j="490pko2004", nsim=100)

## Creating Plot objects
drc.low.2004 <- data.frame(bik.low.drc.sim.2004[0:50])
drc.low.data.2004 <- drc.low.2004$obs - drc.low.2004$theo

drc.med.2004 <- data.frame(bik.med.drc.sim.2004[0:50])
drc.med.data.2004 <- drc.med.2004$obs - drc.med.2004$theo

drc.high.2004 <- data.frame(bik.high.drc.sim.2004[0:50])
drc.high.data.2004 <- drc.high.2004$obs - drc.high.2004$theo

plot(drc.low.data.2004, type="s", main="Clustering of Battle Deaths and Peacekeeping Operations, DRC 2004", ylab="Observed - Expected", xlab="Distance", lty=1, col="black", ylim=c(-1, 10))
lines(drc.med.data.2004, type="s", lty=1, col="blue")
lines(drc.high.data.2004, type="s", lty=1, col="red")
abline(h=0, lty=3, col="gray21")
legend("topleft", inset=0, c("Low BD","Med BD","High BD"), fill=c("black", "blue", "red"), cex=0.75, horiz=TRUE)





_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:R-sig-Geo at r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo