R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Efficient spatial polygon search using kd-trees.

RapidPolygonLookup

R Documentation

Efficient spatial polygon search using kd-trees.

Description

Given spatial partitions such as census blocks, ZIP codes or police district boundaries, we are frequently faced with the need to spatially aggregate data. Unless efficient data structures are used, this can be a daunting task. The operation point.in.polygon() from the package sp is computationally expensive. Here, we exploit kd-trees as efficient nearest neighbor search algorithm to dramatically reduce the effective number of polygons being searched. Points that are left unmapped are put through a linear search to find the associated polygon.

Usage

RapidPolygonLookup(XY, polygons, poly.list = NULL, k = 10, N = nrow(XY), 
    poly.id = "fips", poly.id.colname = "census.block", keep.data = TRUE, 
    verbose = 0)

Arguments

`XY`	data frame containing X-Y or (lon-lat, long-lat, longitude-latitude) columns
`polygons`	polygons to crop and add poly centres
`poly.list`	polygon list with three elements: data, polys, and poly.centers as output from function CropSpatialPolygonsDataFrame()
`k`	maximum number of near neighbours to compute. The default value is set to 10
`N`	number of rows of XY to search
`poly.id`	column name in 'poly.list$data' containing the polygon identifier
`poly.id.colname`	desired column name in the output data frame containing the polygon identifier
`keep.data`	retain polygon list and centers for future referece
`verbose`	level of verbosity

Value

The original points augmented with polygon ID are returned along with the poly centers and other call information

Author(s)

Markus Loecher <markus.loecher@gmail.com> and Madhav Kumar <madhavkumar2005@gmail.com>

Examples

data(sf.crime.2012, envir = environment())
data(sf.polys, envir = environment())
cat(nrow(sf.crime.2012), "rows in SF crime \n")

XY.kdtree <- RapidPolygonLookup(sf.crime.2012[,c("X","Y")], poly.list= sf.polys, 
                                  k= 10, N= 1000, 
                                  poly.id= "fips", poly.id.colname= "census.block", 
                                  keep.data= TRUE, verbose= TRUE)

XY.kdtree.DF <- XY.kdtree$XY
table(XY.kdtree.DF$rank, useNA= "always")
hist(XY.kdtree.DF$rank, xlab = "rank of neighbor")

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(RapidPolygonLookup)
Loading required package: sp
Loading required package: RANN
Loading required package: PBSmapping

-----------------------------------------------------------
PBS Mapping 2.69.76 -- Copyright (C) 2003-2016 Fisheries and Oceans Canada

PBS Mapping comes with ABSOLUTELY NO WARRANTY;
for details see the file COPYING.
This is free software, and you are welcome to redistribute
it under certain conditions, as outlined in the above file.

A complete user guide 'PBSmapping-UG.pdf' is located at 
/home/ddbj/local/lib64/R/library/PBSmapping/doc/PBSmapping-UG.pdf

Packaged on 2015-04-23
Pacific Biological Station, Nanaimo

All available PBS packages can be found at
http://code.google.com/p/pbs-software/

To see demos, type '.PBSfigs()'.
-----------------------------------------------------------


Loading required package: RgoogleMaps
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/RapidPolygonLookup/RapidPolygonLookup.Rd_%03d_medium.png", width=480, height=480)
> ### Name: RapidPolygonLookup
> ### Title: Efficient spatial polygon search using kd-trees.
> ### Aliases: RapidPolygonLookup
> 
> ### ** Examples
> 
> data(sf.crime.2012, envir = environment())
> data(sf.polys, envir = environment())
> cat(nrow(sf.crime.2012), "rows in SF crime \n")
20000 rows in SF crime 
> 
> XY.kdtree <- RapidPolygonLookup(sf.crime.2012[,c("X","Y")], poly.list= sf.polys, 
+                                   k= 10, N= 1000, 
+                                   poly.id= "fips", poly.id.colname= "census.block", 
+                                   keep.data= TRUE, verbose= TRUE)
Mapping points missed in initial run using range-search 
[1] "adding missing polygon range information ..."
> 
> XY.kdtree.DF <- XY.kdtree$XY
> table(XY.kdtree.DF$rank, useNA= "always")

   1    2    3    4    5    6    7    8    9   10 <NA> 
 701  172   58   20   20    9    7    1    8    4    0 
> hist(XY.kdtree.DF$rank, xlab = "rank of neighbor")
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>