Exercise - describing point data

In this exercise we'll create kernel smoothed maps of point location data. This technique produces a smooth raster image of point intensity that is easier to interpret than a map showing point locations alone. By the end of this exercise you should be able to compute a kernel smoothed surface showing the density of cattle holdings in Cheshire and to display that surface in Quantum GIS.

The data for this exercise can be found in the CHESHIRE subfolder of the DATA folder on the course CD-ROM. Copy the contents of the CHESHIRE folder to the working folder on your computer.

File name	What is it needed for?
UK_Cheshire_boundary-BNG.SHP, .SHX, .DBF	ESRI shape file set of Cheshire county boundary.
UK_Cheshire_parish-BNG.SHP, .SHX, .DBF	ESRI shape file set of Cheshire parish boundaries.
UK_Cheshire_FMD-BNG.SHP, .SHX, .DBF	ESRI shape file set with location and FMD details for Cheshire cattle farms.

Table 1: Files required for the describing point data exercise.

This exercise assumes that you have completed the exercise titled Visualising Point Data. Start Quantum GIS, load UK_Cheshire_parish-BNG and UK_Cheshire_FMD-BNG as two layers and set the map units for the study area. Click SETTINGS -OPTIONS. Under the Map Tools tab, set the Preferred distance units to `Kilometers' and the Preferred area units to `Square kilometres'. Click `OK'.

Kernel smoothing

The general approach here is that we will use R to do the kernel smoothing, then take the results from R and load them into Quantum GIS. Although facilities exist within Quantum GIS for kernel smoothing, they lack the flexibility required for good epidemiological analysis.

setwd("D:\\TEMP")
library(rgdal); library(spatstat); library(spatstat.utils); library(maptools)

chpol.shp <- readOGR(dsn = getwd(), layer = "UK_Cheshire_boundary-BNG")
chfmd.shp <- readOGR(dsn = getwd(), layer = "UK_Cheshire_FMD-BNG")

Set the number of grid cells, the bandwidth (1500 metres) and create an observation window:

sigma <- 1500; dimyx <- c(100, 100)
chpol.ow <- as(as(chpol.shp, "SpatialPolygons"), "owin")

Create a ppp object:

pop.ppp <- ppp(x = chfmd.shp$XCOORD, y = chfmd.shp$YCOORD, window = chpol.ow, check = FALSE)

Compute density surface. The density.ppp function returns an estimate of the intensity of the underlying point process, that is the expected number of points per unit area (in this case, the number of points per metre). Express results as the number of points per square kilometre (1 square metre = 0.000001 square kilometres):

pop.den <- density(pop.ppp, sigma = sigma, dimyx = dimyx)
pop.den$v <- pop.den$v / 0.000001

Summarise the holding density estimates:

summary(as.vector(pop.den$v))

The density of holdings ranges from 0.3 to 5.3 per square kilometre. Convert to a SpatialGridDataFrame and set the projection. Plot the density surface:

pop.sgdf <- as.SpatialGridDataFrame.im(pop.den)
proj4string(pop.sgdf) <- CRS("+init=epsg:27700")
spplot(pop.sgdf)

Export the raster image as a GeoTiff file:

writeGDAL(dataset = pop.sgdf, fname = "UK_Cheshire_farm_den.tif", drivername = "GTiff", type = "Float32")

Return to Quantum GIS and add the UK_Cheshire_farm_den.tif as a layer to your project by clicking LAYER - ADD RASTER LAYER. Quantum GIS will add the raster layer as a uniformly coloured grid so some work will be required to make it more interpretable. Right click on the name of the layer and select Properties. Under the Style tab, select `Rendertype' and select singleband pseudocolor. Select Linear as the Color interpolation method and then (at the bottom of the form) set 6 equal interval entries and click Classify. Edit the cutpoints, remembering that holding density in the study area ranged from 0.3 to 5.2 per square kilometre (so 0 to 5 would be a sensible range). Set the colour for a holding density of 0 to 1 to transparent. The Layer Properties dialogue box should look something like Figure 1.

Figure 1: The layer properties form in Quantum GIS showing the Colormap setup for holding density raster surface.

Click Apply then OK. The Properties forms will close and you'll return to your project. Drag the UK_CHESHIRE_parish layer so it is positioned on top of the holding density raster surface. Set it's fill colour to transparent, as shown in Figure 2.

Figure 2: Kernel smoothed surface of farm holding density in the county of Cheshire. The highest concentration of farm holdings is in the south east corner of the county.

Questions

1. Create another surface using a different bandwidth. How does altering the bandwidth change the resulting kernel smoothed surface?
2. Produce a map showing the areas of Cheshire where more than 20 farms per 100 farms per square kilometre were FMD positive.

Answer

Space and time

Kernel smoothed surfaces for case holdings can be constructed and displayed to map users as a contour map (rather than as an image). In Figure 4 contour plots showing the distribution of FMD-positive holdings during each week of the epidemic have been combined into an animated GIF file. This allows the temporal progression of the epidemic to be better appreciated.

Plots of case holdings per unit area lets you appreciate the extent of disease throughout a study area. Plots of case holdings per 100 holdings per unit area lets you distinguish areas where the incidence of disease is high versus areas where the incidence of disease is low.

Figure 4: Point map showing location of FMD positive holdings in the county of Cheshire. The superimposed contour line delineates areas where the intensity of FMD-affected holdings is greater than 10 per 100 holdings per square kilometre. The above series shows that over a period of three weeks infection progresses from the centre of the county to the north eastern corner.

Animated GIF files