Spatio-Temporal Point Pattern Analysis in ArcGIS with R

This post would probably be the last in my series about merging R and ArcGIS. In August unfortunately I would have to work for real and I will not have time to play with R-Bridge any more.
In this post I would like to present a toolbox to perform some introductory point pattern analysis in R through ArcGIS. Basically, I developed a toolbox to perform the tests I presented in my previous post about point pattern analysis. In there, you can find some theoretical concepts that you need to know to understand what this toolbox can do.
I will start by introducing the sample dataset we are going to use, and then simply show the packages available.

Dataset
For presenting this toolbox I am using the same dataset I used for my previous post, namely the open crime data from the UK. For this post I downloaded crimes in the London area from the whole 2015. As you can see from the image below we are talking about more than 950'000 crimes of several categories, all across London.

I also included a polygon shapefile with the area around London and all its boroughs, this should be visible as blue lines around the city. I included this because point pattern analysis requires the user to set the border of the study area, as I mentioned in my previous post.


Spatio-Temporal Subset
The first package I would like to present is a simple spatio-temporal subsetting tool. This is completely based on R but it is basically a more flexible version of the selection tools available in ArcGIS.

Here users can select points based on various parameters at once. For example, they can subset the polygon shapefile, for example here I'm extracting the borough of Ealing, and extract points only for this area. Then they can subset by time, with the same strings I presented in my previous post about a toolbox for time series analysis. Optionally, they can also subset the dataset itself based on some categories. In this example I'm extracting only the drug related crimes, committed in Ealing in May 2015.
It is important to point out that in this first version of the toolbox users can only select one element in the SQL statements. For example here I have "name" = 'Ealing'. In ArcGIS users could also put an AND and then specify another option. However, in the R code I did not put a way to deal with multiple inputs and conditions (e.g. AND, OR) and therefore only one option can be specified.
The result is a new shapefile, plotted directly on the ArcGIS console with the required subset of the data, as shown below:

Spatio-Temporal Centroid
As you may already know, ArcGIS provides a function to calculate the centroid of a point pattern. However, if we wanted to test for changes in the centroid location with time we would need to first subset our data and then compute the centroid. What I did in this package is merge these two actions into one. This package, presented in the image below, loops through the dataset, subsetting the point pattern by time (users can choose between daily, monthly and yearly subsets) and then calculates the centroid for each time unit. Moreover, I also added an option to select the statistics to use between mean, median and mode.

The results for the three statistics are presented below:

Spatio-Temporal Density
This tool calculates the point density for specific regions and time frames by subsetting your dataset. This is something that you may be able to obtain directly from ArcGIS, but users would need to first subset their data and then perform the density analysis, this tool groups those two things into one. Moreover, the package spatstat, which is used in R for point pattern analysis has some clear advantages compared to the tool available in ArcGIS. For example, as I mentioned in my post it provides ways to calculate the best bandwidth for the density estimation. In the script this is achieved using the function bw.ppl, but this can be changed if you need a different method, you just need to replace this function with another. Moreover, as pointed out in this tutorial, ArcGIS does not correct for edge effects.
Working with this package is very similar to the others I presented before:

Users need to specify the input point pattern, then a polygon shapefile for the study area, which can be subset to reduce the area under investigation. Then users can include a temporal subsetting (here I used the string "2015-10/" which means from October to the end of the year, please refer to this post for more info) and subset their data extracting a certain category of crimes. Again here the SQL statements cannot include more than one category.
Finally, users need to provide a raster dataset for saving the density result. This needs to be a .tif file, otherwise in my tests the result did not appear on screen. The output of this script is the image below, for the borough of Bromley and only for robberies:

Spatio-Temporal Randomness
This is another tool to perform a test for spatial randomness, the G function I explained in my previous post, but on a subset of the main dataset. In fact, this test is available in ArcGIS under "Multi-Distance Spatial Cluster Analysis (Ripleys K Function)", but in this case we are again performing it on a particular subset of our data.
The GUI is very similar to the other I presented before:

The only difference is that here users also need to provide an output folder, where the plot created by R will be saved in jpeg at 300 dpi. Moreover, this tool also provides users with the point shapefile created by subsetting the main dataset.
The output for the borough of Tower Hamlets and only for drug related crimes in March 2015 is the plot below:

Spatio-Temporal Correlogram
As the name suggests I develop this tool to calculate and plot a correlogram on a spatio-temporal subset of my data. For this example I could not use the crime dataset, since I do not have a continuous variable in it. Therefore I loaded the dataset of ozone measurements from sensors installed on trams here in Zurich that I used for my post about spatio-temporal kriging. This tool uses the function correlog from the package xts to calculate the correlogram. This function takes several arguments among which an increment, the number of permutations and a TRUE/FALSE flag if data are unprojected or not. These are all data that users will need to input once they use the tool and are additional options in the GUI, which for the other points is more or less identical to what I presented before, except for the selection of the variable of interest:

The result is the image below, which is again saved in jpeg at 300 dpi. As for the spatio-temporal randomness tool, a shapefile with the spatio-temporal subset used to calculate the correlogram is also saved and opened in ArcGIS directly.

Download 
The tool is available, along with the sample data, from my GitHub archive:
https://github.com/fveronesi/Spatio-Temporal-Point-Pattern-Analysis

3D Anisotropy parameters

Hello everyone,

I’m trying to create a method for calculating the 5 parameters for the 3D anisotropy in an automatic way.

Basically I created a loop for analysing the range of the horizontal variogram in every direction and extrapolating the maximum (for the angle p). Then it computes the vertical variogram in the direction of maximum continuity and extrapolates the maximum range (for the angle q).

Subsequently, it computes the vertical variogram for the direction of minimum continuity (angle of max continuity + 90), and extrapolates the maximum (for the angle r).

Ultimately, it calculates the ratio between the range of the maximum continuity and the other 2.

I would like to know what the community think of this idea.

Could you please tell me if it makes sense and, if it does, how to improve it?

So far I wrote the following lines:

library(gstat, pos=4)

library(sp, pos=4)

library(lattice)

c1<-read.asciigrid("dtm.asc")

##a covariate that will be used to specify the cutoff of the variogram

data<-read.table("data.txt",sep="",header=T)

##Data format:

## Lat = latitude

## Lon = longitude

## depth = depth

## value = data value

coordinates(data)=~Lat+Lon+depth

##Anisotropy

##First Loop for angle p

mod<-vgm(var(data$value),"Sph",sqrt(areaSpatialGrid(c1)),0)

write.table(matrix(c("Angle","Range"),ncol=2),"hor.txt",col.names=F,row.names=F)

for(i in 0:180){

angl<-as.data.frame(fit.variogram(variogram(value~1,data,cutoff=sqrt(areaSpatialGrid(c1)),width=sqrt(areaSpatialGrid(c1))/15,alpha=c(i),tol.hor=10),mod,fit.method=6))[2,3]

write.table(matrix(c(i,angl),ncol=2),"hor.txt",append=T,row.names=F,col.names=F)}

hor<-read.table("hor.txt",sep="",header=T)

xyplot(Range~Angle,data=hor,type="l",xlim=0:180,scales=list(x=list(tick.number=15),tck=0.5))

max_range<-subset(hor,hor$Range==max(hor$Range))

##Second Loop for angle q

mod_depth<-vgm(var(data$value),"Mat",max((data@coords)[,3]),0)

write.table(matrix(c("Angle","Range"),ncol=2),"ver1.txt",col.names=F,row.names=F)

for(i in 0:180){

angl<-as.data.frame(fit.variogram(variogram(value~1,data,cutoff=max((data@coords)[,3]),width=max((data@coords)[,3])/length((data@coords)[,3]),alpha=c(max(max_range$Angle)),beta=c(i),tol.hor=10,tol.ver=10),mod_depth,fit.method=6))[2,3]

write.table(matrix(c(i,angl),ncol=2),"ver1.txt",append=T,row.names=F,col.names=F)}

ver1<-read.table("ver1.txt",sep="",header=T)

xyplot(Range~Angle,data=ver1,type="l",xlim=0:180,scales=list(x=list(tick.number=15),tck=0.5))

max_range1<-subset(ver1,ver1$Range==max(ver1$Range))

##Third Loop for angle r

mod_depth<-vgm(var(data$value),"Mat",max((data@coords)[,3]),0)

write.table(matrix(c("Angle","Range"),ncol=2),"ver2.txt",col.names=F,row.names=F)

for(i in 0:180){

angl<-as.data.frame(fit.variogram(variogram(value~1,data,cutoff=max((data@coords)[,3]),width=max((data@coords)[,3])/length((data@coords)[,3]),alpha=c(max(max_range$Angle)+90),beta=c(i),tol.hor=10,tol.ver=10),mod_depth,fit.method=6))[2,3]

write.table(matrix(c(i,angl),ncol=2),"ver2.txt",append=T,row.names=F,col.names=F)}

ver2<-read.table("ver2.txt",sep="",header=T)

xyplot(Range~Angle,data=ver2,type="l",xlim=0:180,scales=list(x=list(tick.number=15),tck=0.5))

max_range2<-subset(ver2,ver2$Range==max(ver2$Range))

##Parameters

p=max(max_range$Angle)

q=max(max_range1$Angle)

r=max(max_range2$Angle)

s=max(max_range$Range)/max(max_range1$Range)

t=max(max_range$Range)/max(max_range2$Range)

R Video Tutorial for Spatial Statistics

Welcome everybody, this is a video tutorial that will try to teach you how to use R for spatial statistics and interpolation.

I’m a PhD student in soil science and in particular my project is about digital soil mapping. During my work I often use R particularly for geostatstical interpolation. Despite R being a very powerful statistical language which has became a benchmark for spatial statistics, it is also a relatively difficult language to learn with one of the steepest learning curve among programming languages.

For this reason I though about sharing what I’ve learned during my PhD with all the R community. I prepared this video tutorial because I think that the easy way of learning R is by examples.

The tutorial is structured as follow: the first two lesson are about the basics of R, how to handle R objects, plotting and saving your work. This part is intended for a beginner user who want to learn the language. Starting from the third lesson I will show examples on how to handle and plotting spatial data and rasters, and in the last two lesson I will show you how to perform an ordinary and a universal kriging in R.

In this tutorial I will show you some example that will hopefully help you learning R in a quicker and easier way. However, this course is not a complete R course, because R has lots and lots of different functions for every branch of science. For this reason at the beginning of the course and at the end of the lesson, if appropriate, I will give you some references if you want to deepen your knowledge.

The link to the course is here: http://www.fabioveronesi.net/rtutorial.html

If you have any suggestion at all on how to improve the course or any other feedback, send an e-mail to: r-course@fabioveronesi.net