Hands-on Exercise 5

In this hands-on experience, I will learn how to use spNetwork package to derive network constrained kernel density estimation (NetKDE), and perform network G-function and k-function analysis.

1. Background Information

In this hands-on experience, I will learn more about the functions of spNetwork package:

• to derive network constrained kernel density estimation (NetKDE), and
• to perform network G-function and k-function analysis

2. Dataset used

• Punggol_St, a line features geospatial data which store the road network within Punggol Planning Area.
• Punggol_CC, a point feature geospatial data which store the location of childcare centres within Punggol Planning Area.

3. Install and load R packages

packages = c('sp', 'rgdal', 'spNetwork', 'tmap')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}

More on the packages used:

• spNetwork provides functions to:
  • perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network
  • build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances
• rgdal provides:
  • bindings to the ‘Geospatial’ Data Abstraction Library (GDAL) (>= 1.11.4)
  • access to projection/transformation operations from the PROJ library.
    
  • In this exercise, rgdal will be used to import geospatial data in R and store as sp objects.
• sp provides:
  • classes and methods for dealing with spatial data in R.
    
  • In this exercise, it will be used to manage SpatialPointsDataFrame and SpatialLinesDataFrame, & for performing projection transformation.
• tmap: functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API

4. Import data

• network is line feature
• childcare is point feature
• Read OGR vector maps into Spatial objects using readOGR of rgdal package

network <- readOGR(dsn="data/geospatial", layer="Punggol_St", verbose = FALSE)
childcare <- readOGR(dsn="data/geospatial", layer="Punggol_CC", verbose = FALSE)

5. Print content of SpatialLineDataFrame & SpatialPointsDataFrame by using str function

str(network)
str(childcare)

From the above results, we can see that the Projected CRS is SVY21. Hence, we can assign the EPSG code using the following code chunk:

6. Assign CRS to OGR data

network <- spTransform(network,CRS("+init=epsg:3414"))
childcare <-spTransform(childcare, CRS("+init=epsg:3414"))

7. Visualising the Geospatial Data

• Good practice to visualise geospatial data after assigning CRS

7.1 Visualise using Base R

plot(network)
plot(childcare,add=T,col='red',pch = 19)

7.2 Visalise using mapping function of tmap package

• Helps us to visualise the geospatial data with high cartographic quality and interactive manner

tmap_mode('view')
tm_shape(childcare)+
  tm_dots() +
tm_shape(network)+
  tm_lines()

tmap_mode('plot')

8. Network Constrained KDE (NetKDE) Analysis using spNetwork package

8.1 Prepare the lixels objects using lixelize_lines of spNetwork

• More on lixelize_lines function: To cut a SpatialLines object into lixels with a specified minimal distance
• Set length of a lixel, lx_length to 700m,
• Set minimum length of a lixel, mindist to 350m.

Note: There is another function called lixelize_lines.mc() which provides multicore support.

lixels <- lixelize_lines(network,700,mindist = 350)
lixels

class       : SpatialLinesDataFrame 
features    : 2645 
extent      : 34038.56, 38882.85, 40941.11, 44801.27  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
variables   : 2
names       : oid.LINK_ID,   oid.ST_NAME 
min values  :  1032897752, EDGEDALE PLNS 
max values  :   998635154,      TRACK 24 

Other important notes:

• After cut, if length of final lixel is shorter than minimum distance, then it is added to previous lixel.
• If NULL, then mindist = maxdist/10.
  
• Segments that are already shorter than the minimum distance are not modified

8.2 Generate line centre points using lines_center() of spNetwork

• Generate a SpatialPointsDataFrame with line center points.
• Points are located at center of line based on the length of the line.

samples <- lines_center(lixels)
samples

class       : SpatialPointsDataFrame 
features    : 2645 
extent      : 34041.34, 38840.31, 41109.22, 44795.36  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
variables   : 3
names       : oid.LINK_ID,   oid.ST_NAME,  Ind 
min values  :  1032897752, EDGEDALE PLNS,    1 
max values  :   998635154,      TRACK 24, 2645 

8.3 Performing NetKDE

densities <- nkde(network, 
                  events = childcare,
                  w = rep(1,nrow(childcare)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

More on the arguments:

• kernel_name indicates quartic kernel is used.
  • Why do we use quartic kernel instead of the other kernels?
    • x
  • Other possible kernel methods supported by spNetwork are: triangle, gaussian, scaled gaussian, tricube, cosine, triweight, epanechnikov or uniform.
• method indicates simple method is used to calculate the NKDE. More on the possible methods supported by spNetwork:
  • simple: the distances between events and sampling points are replaced by network distances, and the formula of the kernel is adapted to calculate the density over a linear unit instead of an areal unit.
  • discontinuous: Equally “divides” the mass density of an event at intersections of lixels.
  • continuous: Divide the mass of the density at intersection but adjusts the density before the intersection to make the function continuous

8.4 Visualising NetKDE

8.4.1 Insert computed density values into samples & lixels objects

samples$density <- densities
lixels$density <- densities

8.4.2 Rescale density values

• Why we have to rescale? : Since svy21 projection system is in meter, the computed density values are very small i.e. 0.0000005.

# rescaling to help the mapping
samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

8.4.3 Plot map after rescaling

tmap_mode('view')
tm_shape(lixels)+
  tm_lines(col="density")+
tm_shape(childcare)+
  tm_dots()

tmap_mode('plot')

What can we see from map above:

• Road segments with relatively higher density of childcare centres (darker color)
• Road segments with relatively lower density of childcare centres (lighter color)

9. Network Constrained G- and K-Function Analysis

• Here, we will perform complete spatial randomness (CSR) test by using kfunctions() of spNetwork package.
  
• The null hypothesis is defined as: Ho: The observed spatial point events (i.e distribution of childcare centres) are UNIFORMLY DISTRIBUTED over a street network in Punggol Planning Area.
• The CSR test is based on the assumption of the binomial point process which implies the hypothesis that the childcare centres are randomly and independently distributed over the street network.
• If this hypothesis is rejected, we may infer that the distribution of childcare centres are spatially interacting and dependent on each other; as a result, they may form nonrandom patterns.

kfun_childcare <- kfunctions(network, 
                             childcare,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)

The output of kfunctions() is a list with the following values:

• plotkA, a ggplot2 object representing the values of the k-function
• plotgA, a ggplot2 object representing the values of the g-function
• valuesA, a DataFrame with the values used to build the plots

Just for my understanding:

9.1 kfunctions g-function

kfun_childcare$plotg

9.2 kfunctions values

kfun_childcare$values

         obs_k     lower_k    upper_k      obs_g    lower_g   upper_g
1      0.00000     0.00000     0.0000    0.00000    0.00000  200.7234
2      0.00000     0.00000   417.7217   72.33276   16.27487  417.7217
3     72.33276    72.33276   707.0527    0.00000   72.33276  506.3293
4     72.33276   216.99828  1213.3820   72.33276   32.54974  690.7779
5    361.66380   522.60419  1614.8288  289.33104   72.33276  707.0527
6    650.99483   884.26798  2226.0407  216.99828  160.94039  980.1089
7    795.66035  1301.98967  3173.5998  289.33104  216.99828 1141.0493
8   1229.65691  1712.47807  4018.0848  578.66207  305.60591 1213.3820
9   2097.65002  2169.98278  5119.3510 1084.99139  305.60591 1430.3803
10  3182.64140  2837.25248  6188.0675 1084.99139  377.93867 1446.6552
11  4195.30003  3464.73917  7658.2309  650.99483  522.60419 1864.3769
12  4773.96211  4164.55861  9434.0001 1301.98967  594.93694 2282.0986
13  6654.61385  5168.17564 11627.4910 2097.65002  795.66035 2009.0424
14  8390.60007  6108.50151 13365.2856 1663.65346  940.32587 2482.8220
15 10054.25353  7482.82394 15723.3335 1591.32070  916.81772 2788.4279
16 12079.57079  8535.26559 17797.4754 1808.31898 1012.65863 2660.0372
17 14032.55528  9974.68749 19582.2862 2097.65002 1084.99139 2716.0951
18 16202.53806 11300.18531 22057.8749 2603.97933  956.60074 3166.3665
19 18661.85187 12547.92540 24676.3208 2459.31381 1374.32242 2989.1513
20 21338.16396 14386.98581 27249.5587 2676.31209 1334.53941 2965.6431
21 24014.47605 15786.62470 29218.8181 2712.47847 1021.70022 3238.6993
   distances
1          0
2         50
3        100
4        150
5        200
6        250
7        300
8        350
9        400
10       450
11       500
12       550
13       600
14       650
15       700
16       750
17       800
18       850
19       900
20       950
21      1000

Other notes about the arguments used:

• lines: A SpatialLinesDataFrame with the sampling points. The geometries must be a SpatialLinesDataFrame (may crash if some geometries are invalid).
• points: A SpatialPointsDataFrame representing the points on the network. These points will be snapped on the network.
• start: A double, the start value for evaluating the k and g functions.
• end: A double, the last value for evaluating the k and g functions.
• step: A double, the jump between two evaluations of the k and g function.
• width: The width of each donut for the g-function.
• nsim: An integer indicating the number of Monte Carlo simulations required. In the above example, 50 simulation was performed. Note: most of the time, more simulations are required for inference
• resolution: When simulating random points on the network, selecting a resolution will reduce greatly the calculation time. When resolution is null the random points can occur everywhere on the graph. If a value is specified, the edges are split according to this value and the random points are selected vertices on the new network.
• conf_int: A double indicating the width confidence interval (default = 0.05).

10. Visualise ggplot2 object of k-function

kfun_childcare$plotk

From the plot above:

• blue line is the empirical network K-function of the childcare centres in Punggol planning area
• gray envelop represents the results of the 50 simulations in the interval 2.5% - 97.5%.
• Because blue line between the distance of 250m-400m are below the gray area, we can infer that the childcare centres in Punggol planning area resemble regular pattern at the distance of 250m-400m.