Author: Serge Rey sjsrey@gmail.com
Windows play several important roles in the analysis of planar point patterns. As we saw in the introductory notebook, the area of the window can be used to develop estimates of the intensity of the point pattern. A window also defines the domain for the point pattern and can support corrections for so-called edge effects in the statistical analysis of point patterns. However, there are different ways to define a window for a point pattern.
This notebook provides an overview of how to work with windows and covers the following:
We will first continue on with an example from the introductory notebook. Recall this uses 200 randomly distributed points within the counties of Virginia. Coordinates are for UTM zone 17 N.
In [1]:
import pysal.lib as ps
import numpy as np
from pysal.explore.pointpats import PointPattern
In [2]:
f = ps.examples.get_path('vautm17n_points.shp')
fo = ps.io.open(f)
pp_va = PointPattern(np.asarray([pnt for pnt in fo]))
fo.close()
pp_va.summary()
From the summary method we see that the Bounding Rectangle is reported along with the Area of the window for the point pattern. Two things to note here.
First, the only argument we passed in to the PointPatterns constructor was the array of coordinates for the 200 points. In this case PySAL finds the minimum bounding box for the point pattern and uses this as the window.
The second thing to note is that the area of the window in this case is simply the area of the bounding rectangle. Because we are using projected coordinates (UTM) the unit of measure for the area is in square meters.
The window is an attribute of the PointPattern. It is also an object with its own attributes:
In [3]:
pp_va.window.area
Out[3]:
In [4]:
pp_va.window.bbox
Out[4]:
The bounding box is given in left, bottom, right, top ordering.
In [5]:
pp_va.window.centroid
Out[5]:
In [6]:
pp_va.window.parts
Out[6]:
The parts attribute for the window is a list of polygons. In this case the window has only a single part and it is a rectangular polygon with vertices listed clockwise in closed cartographic form.
A window has several basic geometric operations that are heavily used in some of the other modules in the the Point package. Most of this is done under the hood and the user typically doesn't see this. However, there can be times when direct access to these method can be handy. Let's explore.
The window supports basic point containment checks:
In [7]:
pp_va.window.contains_point((623277.82697965798, 4204412.8815969583))
Out[7]:
This also applies to sequences of points:
In [8]:
pnts = ((-623277.82697965798, 4204412.8815969583),
(623277.82697965798, 4204412.8815969583),
(1000.01, 200.9))
In [9]:
pnts_in = pp_va.window.filter_contained(pnts)
pnts_in
Out[9]:
In [10]:
parts = [[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0)],
[(11.,11.), (11.,20.), (20.,20.), (20.,11.)]]
holes = [[(3.0,3.0), (6.0, 3.0), (6.0, 6.0), (3.0, 6.0)]]
We will plot this using matplotlib to get a better understanding of the challenges that this type of window presents for statistical analysis of the associated point pattern.
In [11]:
%matplotlib inline
import matplotlib.pyplot as plt
p0 = np.asarray(parts[0])
plt.plot(p0[:,0], p0[:,1])
plt.xlim(-10,20)
t = plt.ylim(-10,20) # silence the output of ylim
Not, quite what we wanted, as the first part of our multi-part polygon is a ring, but it was not encoded in closed cartographic form:
In [12]:
p0
Out[12]:
We can fix this with a helper function from the window module:
In [13]:
from pysal.explore.pointpats.window import to_ccf
print(parts[0])
print(to_ccf(parts[0])) #get closed ring
In [14]:
from pysal.explore.pointpats.window import to_ccf
p0 = np.asarray(to_ccf(parts[0]))
plt.plot(p0[:,0], p0[:,1])
plt.xlim(-10,20)
t=plt.ylim(-10,20)
Now we can print all the rings composing our window: two exterior rings, and one hole:
In [15]:
for part in parts:
part = np.asarray(to_ccf(part))
plt.plot(part[:,0], part[:,1], 'b')
for hole in holes:
hole = np.asarray(to_ccf(hole))
plt.plot(hole[:,0], hole[:,1], 'r')
plt.xlim(-10,30)
t = plt.ylim(-10,30)
The red hole is associated with the first exterior ring.
With this visual representation, consider the problem of testing whether or not this multi-part window contains one or more points in a sequence:
In [16]:
pnts = [(12,12), (4,4), (2,2), (25,1), (5,20)]
for pnt in pnts:
plt.plot(pnt[0], pnt[1], 'g.')
for part in parts:
part = np.asarray(to_ccf(part))
plt.plot(part[:,0], part[:,1], 'b')
for hole in holes:
hole = np.asarray(to_ccf(hole))
plt.plot(hole[:,0], hole[:,1], 'r')
plt.xlim(-10,30)
t = plt.ylim(-10,30)
Of the five points two are clearly outside of both of the exterior rings. The three remaining points are each contained in one of the bounding boxes for an exterior ring. However, one of these points is also contained in the hole ring, and thus is not contained in the exterior ring associated with that hole.
We can create a Window object from the parts and holes to demonstrate how to evaluate these containment checks.
In [17]:
from pysal.explore.pointpats import Window
window = Window(parts, holes)
In [18]:
window.parts
Out[18]:
In [19]:
window.holes
Out[19]:
In [20]:
window.bbox
Out[20]:
In [21]:
window.area
Out[21]:
In [22]:
pnts = [(12,12), (4,4), (2,2), (25,1), (5,20)]
for pnt in pnts:
plt.plot(pnt[0], pnt[1], 'g.') #plot the five points in green
for part in parts:
part = np.asarray(to_ccf(part))
plt.plot(part[:,0], part[:,1], 'b') #plot "parts" in blue
for hole in holes:
hole = np.asarray(to_ccf(hole))
plt.plot(hole[:,0], hole[:,1], 'r') #plot "hole" in red
from pysal.explore.pointpats.window import poly_from_bbox
poly = np.asarray(poly_from_bbox(window.bbox).vertices)
plt.plot(poly[:,0], poly[:,1], 'm-.') #plot the minimum bounding box in magenta
plt.xlim(-10,30)
t = plt.ylim(-10,30)
Here we have extended the figure to include the bounding box for the multi-part window (in cyan). Now we can call the filter_contained method of the window on the point sequence:
In [23]:
pin = window.filter_contained(pnts)
pin
Out[23]:
This was a lot of code just to illustrate that the methods of a window can be used to identify topological relationships between points and the window's constituent parts. Let's turn to a less contrived example to see this in action.
Here we will make use of PySAL's shapely extension to create a multi-part window from the county shapefile for Virgina.
In [24]:
from pysal.lib.cg import shapely_ext
import numpy as np
from pysal.explore.pointpats.window import poly_from_bbox, as_window, Window
import pysal.lib as ps
%matplotlib inline
import matplotlib.pyplot as plt
In [25]:
va = ps.io.open(ps.examples.get_path("vautm17n.shp")) #open "vautm17n" polygon shapefile
polys = [shp for shp in va]
vapnts = ps.io.open(ps.examples.get_path("vautm17n_points.shp")) #open "vautm17n_points" point shapefile
points = [shp for shp in vapnts]
In [26]:
print(len(polys))
The county shapefile vautm17n.shp has 136 shapes of the polygon type. Some of these are composed of multiple-rings and holes to reflect the interesting history of political boundaries in that State.
Fortunately, with our window class we can handle these. We will come back to this shortly.
First we are going to build up a realistic window for our point pattern based on a cascaded union made possible via Shapely through the PySAL shapely extension.
In [27]:
cu = shapely_ext.cascaded_union(polys)
This creates a PySAL Polygon:
In [28]:
type(cu)
Out[28]:
We can construct a Window from this polygon instance using the helper function as_window:
In [29]:
w = as_window(cu)
In [30]:
w.holes
Out[30]:
In [31]:
len(w.parts)
Out[31]:
The window has three parts consisting of the union of mainland counties and two "island" parts associated with Accomack and Northampton counties and has no holes.
Since this a window, we can access its properties:
In [32]:
w.bbox
Out[32]:
In [33]:
w.centroid
Out[33]:
In [34]:
w.contains_point(w.centroid)
Out[34]:
So the centroid for our new window is contained by the window. Such a result is not guaranteed as the geometry of the window could be complex such that the centroid falls outside of the window.
Let's continue on with a more interesting query. Since we know the window centroid is contained in the Window, we can find which individual county contains the centroid.
Our strategy is a simple one to illustrate the useful nature of the Window. We will create a sequence of Windows, one for each county and use them to carry out a containment test.
In [35]:
#create a window for each of the individual counties in the state
windows = [as_window(county) for county in polys]
In [36]:
#check each county for containment of the window's centroid
cent_poly = [ (i, county) for i,county in enumerate(windows) if county.contains_point(w.centroid)]
cent_poly
Out[36]:
In [37]:
i, cent_poly = cent_poly[0]
In [38]:
cent_poly.bbox
Out[38]:
What we did here was create a window for each of the individual counties in the state. With these in hand we checked each one for containment of the window's centroid. The result is we see the window (count) with index 67 is the only one that contains the centroid point.
The point of this exercise is not to use an inefficient brute force exhaustive search to find this county. There are more efficient spatial indices in PySAL that we could use for such a query. Rather, we wanted to explicitly check each window to ensure that only one contained the centroid.
As we will see in elsewhere in this series of notebooks, this type of decomposition can support highly flexible types of spatial analysis.
Returning to the central use of Windows, we saw in the introductory notebook that the area of the Window is used to form the estimate of intensity for the point pattern:
In [39]:
f = ps.examples.get_path('vautm17n_points.shp') #open "vautm17n_points" point shapefile
fo = ps.io.open(f)
pnts = np.asarray([pnt for pnt in fo])
fo.close()
pp_va = PointPattern(pnts)
pp_va.summary()
Here the default is to form the minimum bounding rectangle and use that as the window for the point pattern and, in turn, to implment the intesity estimation.
We can override the default by passing a window object in to the constructor for the point pattner. Here we use our window that was formed from the county cascading union above:
In [40]:
pp_va_union = PointPattern(pnts, window=w)
pp_va_union.summary()
Here, the window is redefined. Thus, window related attributes Area of window and Intensity estimate for window are changed. However, the Bounding rectangle remains unchanged since it is not relavant to the definition of window.
Close examination of the summary report from reveals that while the bounding rectangles for the two point pattern instances are identical (as they should be), the area of the windows are substantially different:
In [41]:
pp_va.window.area / pp_va_union.window.area
Out[41]:
as are the intensity estimates:
In [42]:
pp_va.lambda_window / pp_va_union.lambda_window
Out[42]:
In [ ]: