stampr : Spatial-Temporal Analysis of Moving Polygons in R

The R package stampr implements functions for analyzing movement in mapped polygon data. Methods described in this paper include deriving change events based on spatial relationships, plotting change events, summarizing measures of distance and direction of movement, characterizing changes in polygon shape changes, and characterizing sequences of polygons over time using graphs. Two examples are used to demonstrate the core functionality available in the stampr package.


Introduction
Spatial and spatial-temporal analytical methods in many fields have seen considerable development in recent years. However most methods capture only spatial correlation or covariance structures over polygonal lattice data, or spatial characteristics of point patterns. For polygon and line spatial data, there are fewer methods available for spatial-temporal analysis. In this paper we consider the analysis of spatial-temporal polygon data (termed moving polygons; see Robertson, Nelson, Boots, and Wulder 2007 for example). The analysis of moving polygons is widely applicable to the environmental sciences, for example in the spread of zoonotic infectious disease (Morris, Blackburn, Talibzade, Kracalik, Ismaylova, and Abdullahyev 2013), wildlife habitat analysis (Smulders, Nelson, Jelinski, Nielsen, Stenhouse, and Laberee 2012), and spatial changes in land use over time (Mizutani 2012). In geographic information systems (GIS), natural phenomena tend to be geometrically represented as either vector points, lines, polygons (i.e., the object) model, or raster surfaces (i.e., thefield model). Yet many phenomena do not readily fall into these categories, such as wildfires, hurricanes, or ecological hotspots, all of which exhibit both object and field characteristics. Analytical methods for characterizing space-time change in such phenomena have been poorly developed in part due to these issues of spatial representation. For polygon-change analysis, methods and compu-tational tools for examining changes in spatial polygons through time remain limited, most notably in their accessibility to researchers through software tools.
With the growing availability of multi-temporal spatial polygon data, there is a clear need for tools and techniques for space-time analysis of polygons. As such, we developed a new package in the statistical computing environment R (R Core Team 2017), the stampr package (Long and Robertson 2018). Package stampr is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=stampr and provides a unified implementation of old and new approaches to spatial-temporal analysis of moving polygons. The stampr package extends the core STAMP methodology presented by Robertson et al. (2007, hereafter referred to as RNBW) for categorizing space-time change in moving polygons. Within the stampr package, we also include a suite of methods for assessing distance and directional relationships between moving polygons. We extend the STAMP methodology to include polygon shape metrics in order to facilitate more in-depth analysis of polygon movement and change, and introduce a graph representation of space-time change trajectories for examining temporal patterns of spatial change. In all the work here, we consider discrete temporal changes only and leave interpolation of continuous change to future research.
The availability of R-based spatial classes and analysis methods has increased substantially in recent years, owing heavily to the developments of a relatively small number of contributors working on spatial packages in R. For more details on R packages available for handling analyzing spatio-temporal data see the corresponding CRAN Task View (Pebesma 2017). The foundational spatial classes in R are from the sp package (Pebesma and Bivand 2005) for vector data and the raster package (Hijmans 2017) for raster data. The development of the rgeos package (Bivand and Rundel 2017), which implements spatial-overlay operations in R, has been integral to the development of the functionality presented here. For many users who previously relied on proprietary software, the R statistical computing environment has become the preferential free and open-source option as it now possesses functionality, historically available only in commercial GIS packages. R is advantageous for many reasons, but notably due to its wide-reaching user-base and the ability to perform statistical analysis with spatial data through their use of the 'data.frame' class. With this in mind, the R software environment represents an ideal computing environment for reaching potential users of STAMP in both the natural and social sciences.

STAMP framework for space-time analysis of polygons
Metric spatial relationships are not well defined between multiple polygons (Liu and Deng 2002). As such, polygon spatial relations and resultant spatial analyses are typically confined to topological relationships defined by point-set topology (boundary and interior areas of intersection and non-intersection; see Egenhofer and Franzosa 1991). The complexity of polygon features makes characterizing the spatial relationships between polygons a challenging analytical task. Furthermore, the lack of a more sophisticated framework for the spatial analyses of polygonal data and the paucity of analytical tools in both commercial and open source GIS software limits the full exploitation of the geometric properties of the data. In the context of moving polygons, metric relations (e.g., Hausdorff distance), are useful to describe changes in edges or polygon centroids. When aggregated over multiple time periods, changes in distance and direction of movement are descriptive of the change over time in the underlying spatial process.
The two basic metric spatial relations between moving polygons are the distance (related to velocity) and direction of movement. Distance and direction may be computed between polygon centroids or edges. A confounding issue in quantifying distance and directional relationships in moving polygons is that changes to polygon shape and size occur simultaneously with movement (e.g., a wildfire or hurricane). Interestingly, the application domains where analysis of moving polygons is typically required often lie at the interface of discrete objects and continuous fields, described as field-objects (Cova and Goodchild 2002). Field-objects can be dynamic in location, shape, orientation, internal homogeneity, and size, making evaluation of metric relations between them cumbersome due to interacting simultaneous changes in these properties.

Spatial-temporal analysis of moving polygons (STAMP)
Building upon the original works of Sadahiro and Umemura (2001), RNBW laid the groundwork for a more quantitative approach to analyzing space-time change in moving polygons. RNBW implemented an event-based framework for categorizing polygon movement events based on the intersection of t1 and t2 polygons, and a distance threshold (dc) used to define whether two (or more) polygons are related (termed groups) between t1 and t2. Grouped polygons are then evaluated within a hierarchical system used to delineate various polygon movement events in order to semantically describe polygon movements. At the broadest hierarchical level (level 1) STAMP categorizes three event types, stable (areas that exist in both times), generation (areas in t2 only), and disappearance (areas in t1 only). At level 2, categories include stable (as before), expansion (generation events adjacent to stable events), contraction (disappearance events adjacent to stable events), and displaced disappearance and displaced generation events. At level 3, the distance threshold is used to further refine the STAMP displacement designations. Stable, contraction, and expansion events remain unchanged, but disjointed disappearance and generation events are further classified into sub-events based on their proximity and relationship to the other polygon event types (see Table 2 in RNBW for a more detailed description). Finally, level 4 of the framework describes multi-polygon events which include union of two or more polygons into one (i.e., merging), division of one polygon into two or more (i.e., splitting), or both occurring simultaneously. Using a simple two-time period synthetic dataset available in the stampr package, we will demonstrate the variety of scenarios that can be studied under the STAMP event framework in the stampr package.
We begin by running the stamp function in order to generate a change object of class 'SpatialPolygonsDataFrame' that describes change according to the STAMP event framework. These data represent simulated fire perimeter polygons from two time periods. R> library("stampr") R> data("fire1", package = "stampr") R> data("fire2", package = "stampr") The function requires temporally unique identifiers to be stored in a column called ID in the 'SpatialPolygonsDataFrame' object passed to the function. We will create this prior to running the stamp function R> fire1$ID <-1:nrow(fire1) R> fire2$ID <-(max(fire1$ID) + 1):(max(fire1$ID) + nrow(fire2)) R> ch <-stamp(fire1, fire2, dc = 1, direction = FALSE, distance = FALSE, + shape = FALSE) The order of arguments in the stamp function is important; it must be an object of class 'SpatialPolygonsDataFrame' from t1, followed by t2. The returned object ch has class 'SpatialPolygonsDataFrame' with specific fields describing the events and changes between the two time periods. which in Figure 2 shows that one polygon group has higher proportional contraction, while three others have higher proportional expansion. A balance of expansion and contraction would indicate movement across the landscape. We will now investigate further the spatial relationships between polygons from neighboring time periods.

Number, area, and shape indices
RNBW first characterized changes between two polygons using global change indices which can be used to describe the changes in number and area of polygons contained in t1 and t2 polygon groups (see Table 4 in RNBW). Further, we identify a set of shape indices, mostly from the spatial ecology literature (McGarigal and Marks 1995), that are useful for measuring various properties of polygon shapes. Shape indices implemented include perimeter-area ratio, fractal dimension, shape index, and linearity index (information on how each index is calculated can be found in the package documentation, or see McGarigal and Marks 1995). Each index can be used to evaluate shape-related change alongside the distance and directional indices in order to better contextualize the changing dynamics of moving polygons.
Inherent in these event-based change descriptors are more complex changes over time, for example splitting and merging of polygons. Polygons that split or merge through time induce complexities into their distance and directional relationships due to changing the spatial context within which two (or more) polygons are compared. For example, Stell and Worboys (2011) provided recently the theoretical basis for computing relations for splitting and merging polygons using adjacency trees, with focus on the relationship between polygons (in the foreground) and a background region. We will extend the implementation to handle distance and direction between level 4 events in future releases.

Spatial relationships of polygons: Distance and direction
Methods for quantifying distance relationships in polygon sets have developed largely from known analytical problems encountered in computational geometry (e.g., De Berg, Van Kreveld, Overmars, and Schwarzkopf 2008) and other disciplines. In practice, the calculation of distance indices between polygons can be relatively straightforward, and is often simplified by decomposing a polygon (consisting of vertices and edges) into a point-set (i.e., just the vertices) in order to perform calculations. Even more simplistic methods represent the polygon as a single point (e.g., centroid) in order to calculate distances. This facilitates relatively straightforward algebra for computing distance relationships. However, given the complex spatial shapes that can arise in polygon analysis, simple distance metrics can be confounded by changes in polygon shape.
Conversely, algorithms for performing polygon directional analysis need to take advantage of the specific properties associated with polygons and their directional relationships. The formal properties of a directional relationship between two polygons in the plane are identified by Peuquet and Zhang (1987) as follows: I. it is binary; II. it is semantic inverse; III. the area of acceptance of a specified direction increases with distance; IV. the area of acceptance of a specified direction increases with the dimension in the facing direction of the reference polygon in relation to the target polygon.
The first property states that a directional relationship can only involve two polygons, and the second states that the reverse of a given relation is also true. For example, in Figure 3a, if polygon A is north of polygon B, then polygon B is south of polygon A. However, polygons with complex shapes and directional relationships may challenge this convention. For example in Figure 3b polygon A is both north and east of polygon B. The third and fourth properties require more careful consideration for implementation and interpretation.
As stated previously in property 3 (distance property), the distance between two polygons impacts the interpretation of directional relationships between them. For example, in Figure 3b, polygon A is east of polygon B, but the entirety of polygon B is also south of polygon A. This situation is made more complex when different types of distances between two simple polygons are considered.
Polygon shape, size and orientation are additional spatial properties associated with polygons shown to impact distance and directional relationships (Miller and Wentz 2003). These are summarized in property 4 (shape property) which states that the directional relationship is impacted by the ratio of the facing dimensions between the reference and target polygons ( Figure 3c). When this ratio is high, the directional relationship is high, and when this ratio is low, the directional relationship is more ambiguous. This is illustrated in Figure 3c, where the directional relation of polygon A being north of polygon B is stronger than the relation of polygon C being south of polygon B. Figure 3c combines these properties of distance and shape  to demonstrate their combined impact on directional relationships. In stampr we provide two polygon distance methods that facilitate the calculation of distance measurements between moving polygon objects.

Centroid distance
stamp.distance(ch, dir.mode = "Centroid") The centroid distance is the simplest and most straightforward index for quantifying the movement distance of two polygons. Simply put, the centroid distance is the spatial distance between the centroid of two polygons (Figure 4a). The centroid distance can be evaluated between the grouped t1 and t2 polygons, or across the individual STAMP events.

Hausdorff distance
stamp.distance(ch, dir.mode = "Hausdorff") The Hausdorff distance is a metric distance used in a variety of applications for comparing the distance between two point-sets (Figure 4b). In application to measuring polygon distance implemented here, the Hausdorff distance measures the maximum distance between the corresponding edges of two polygons. The Hausdorff distance is more sensitive to changes in shape than the centroid distance. Other applications have also adapted the Hausdorff distance as a metric with varied success, for example in evaluating the similarity of moving point objects (Zhang, Huang, and Tan 2006;Shao, Cai, and Gu 2010). The Hausdorff distance is defined as with i and j the indices along the boundaries of the two polygons P 1 and P 2 respectively, and d a distance operator (e.g., the Euclidean distance). The Hausdorff distance can be interpreted as a measure of separation, the maximum distance separating the two polygons. The Hausdorff distance can be applied to the grouped t1 and t2 polygons, or to individual STAMP events. Variations of the Hausdorff distance can also be used, for example the forward max-min distance from t1 to t2 polygons.
We provide a suite of four methods for describing directional change with moving polygons. Some of these methods were compared over polygons in previous work (Robertson and Nelson 2008).

Centroid angle
stamp.direction(ch, dir.mode = "CentroidAngle") The simplest and most straightforward method for polygon directional analysis is the centroid angle method, which measures the angle between two polygon centroids ( Figure 5a). In the context of STAMP, one can measure the distance between the t1 polygon centroid, and the centroid of each stamp movement event to investigate the direction of each movement event type. Similarly, one can measure the distance between the centroid of grouped t1 and t2 polygons to get a single index of direction for each polygon group.

Cone method
stamp.direction(ch, dir.mode = "ConeModel") The cone method for directional analysis was proposed by Peuquet and Zhang (1987). In this approach, an inverted triangle originates from the centroid of the t1 polygon, typically with an interior angle of 90 • for the ndir = 4 cardinal directions (Figure 5b) or 45 • for the ndir = 8 cardinal directions. The number of cones can be modified to any number, with the interior angle then calculated appropriately. The area of the t2 polygon lying inside each of the cones is computed (and can be converted to percentages if desired). Often it is useful to identify the cone with the maximum area to identify the main directional relation. The cone method clearly satisfies the distance property (property III earlier), as distance from the reference polygon increases, the area of acceptance, or area associated with each possible direction also increases due to the shape of the cone. The cone method is straightforward to compute and interpret, and can be easily visualized. The cone method can be evaluated on the t1 and t2 polygon groups, or on the individual STAMP event polygons.

Modified cone method
stamp.direction(ch, dir.mode = "ModConeModel") The cone method was modified by RNBW to better accommodate irregularly shaped polygons and scenarios where multiple polygons overlap. Like the cone method, the modified cone method starts with the centroid of the t1 polygons (Figure 5c). The bounding box around the union of the t1 and t2 polygons is then used as a reference for extending either ndir = 4 or ndir = 8 uneven cones. If ndir = 4 cones are extended outward to the corners of the bounding box. If ndir = 8 then cones are extended outward to the intersection points along the edge of each side of the bounding box. Thus, while still incorporating the distance property, this approach also attempts to incorporate the shape property by adjusting the geometry of the cone based on the shapes and orientation of the reference and target polygons. RNBW made further modifications using generalized Voronoi diagrams for partitioning special cases, such as when there are multiple stable events within a single STAMP group. The modified cone method can be applied to t2 polygons as a group, or on the individual STAMP event polygons.

Minimum bounding rectangle method
stamp.direction(ch, dir.mode = "MBRModel") The minimum bounding rectangle (MBR) method (Skiadopoulos, Giannoukos, Sarkas, Vassiliadis, Sellis, and Koubarakis 2005) computes the area of movement in 9 quadrants based on the geometry of the t1 polygon. The bounding box of the t1 polygon is used to delineate the central quadrant of an unevenly spaced 3×3 lattice (Figure 5d). Each border of the bounding box is extended outward (ad infinitum) to delineate the 8 outer quadrants for each of the 8 cardinal directions. The area of the t2 polygon in each of the 9 quadrants is calculated. The MBR method can be evaluated on t2 polygons as a group, or on the individual STAMP event polygons.

Graph analysis of polygon sequences
Linking polygons from neighboring time periods through their spatial relationships affords the ability to create connected sequences of polygons that are part of a common change event (or set of events). In Sadahiro (2001), these sequences are termed spatio-temporal histories. Sadahiro (2001) also introduced the graph representation of these sequences through the directed graph, where nodes represented polygons and edges a spatial relationship between polygons at different time periods. The directed graph of polygon change events can be used to summarize changes in size, distance, direction, or shape change over many time periods. As well, we have recently explored applying graph measures as new tools for characterizing different polygon change processes, an approach discussed in Robertson (2016).

Data
Two datasets are included in the stampr package to facilitate testing and understanding of the methods. These are described below.

Hurricane Katrina
The Hurricane Katrina data were obtained from the NOAA H*Wind data product, which provides a gridded measurement of wind speed covering the hurricane region at 3 hour time steps (Powell, Houston, Amat, and Morisseau-Leroy 1998, http://www.rms.com/models/ hwind/legacy-archive). Hurricanes can be categorized using the Saffir-Simpson wind scale, which is based on wind speed. The Saffir-Simpson wind scale ranges from tropical storm (≥ 39 mph) to class 5 hurricane (≥ 157 mph). To delineate polygon boundaries of hurricane Katrina, we identified the 39 mph isotach (line of equal wind speed) from the NOAA H*Wind data, which represents wind-speeds in a regular grid (raster) format. Data are available for every three hours from 21:00 on Friday 2005-08-26 to 21:00 on Monday 2005-08-29. For simplicity, we included the eye of the hurricane as contained within the hurricane polygon, despite the fact that this area of the hurricane will typically have a much lower wind speed. Thus, the hurricane Katrina dataset represents the hurricane boundary polygon at each of 33 time points covering this 3-day period, along with the trajectory of the eye of the hurricane ( Figure 6).

Mountain pine beetle hotspots
The mountain pine beetle dataset were obtained from helicopter GPS mapping of clusters of pine trees infested with mountain pine beetle in western Canada. Mountain pine beetles are a tree-killing insect that emerge from egg galleries in the phloem tissue of the tree's bark each year to attack new trees (Safranyik and Carroll 2006). As the tree begins to die, the tree foliage turns red, providing an indication of where the beetles were the previous year (Wulder, Dymond, White, Leckie, and Carroll 2006). In the late 1990s and early 2000s, the largest mountain pine beetle outbreak ever recorded occurred in western Canada. The point observations of attacked tree-clusters were converted to "hotspot" polygons through kernel density thresholding in order to track how the infestation is spreading from year-toyear (Nelson and Boots 2008). Through spatial-temporal analysis of hotspot movement we can better understand processes driving the epidemic.

Hurricane Katrina analysis
We performed STAMP analysis with the hurricane Katrina data, computing the STAMP events along with each distance index described above. We demonstrate how STAMP analysis can inform on polygon-movement for a spatial process changing over time. With the hurricane Katrina data we can see distinct periods of generation (early Saturday and Sunday) corresponding with the times when the change in polygon area is greatest (steep positive slope of the Katrina polygon area line Figure 7). We also identify that disappearance is largest as hurricane Katrina dissipates late Monday, and the size of the storm area begins to drop.
Comparing between the two polygon-based distance metrics (centroid and Hausdorff distance) we notice that in general the Hausdorff distance is larger than the centroid distance ( Figure 8). Directional analysis reveals that the centroid angle method and the area-based direction methods produce different outputs, leading to differing inferences in some cases. We compare the centroid angle with the area-weighted mean direction method of the three area-based measures for hurricane Katrina at each time-frame (Figure 9).

MPB analysis
The Katrina example shows how multi-temporal polygon analysis can be performed using the stampr package for a single event changing over time. The MPB data provides an example of stampr application for analysis of many polygons, i.e., 711 polygons over 8 time periods. The MPB data represent two infestations, which are distinguished in the REGION variable in the dataset. We have used rose diagrams to identify the dominant directions of expansion and contraction events which correspond to different aspects of movement. We investigated  directional trends in the north and south independently, examining how the directional distribution of expansion and contraction events changed over time (Figures 10 and 11). We see that the dominant direction was in the east-west direction in both regions, likely resulting from wind patterns. We also see a significant southward expansion event in change period 6 in the southern study area (Figure 11).
We visualize the spatial-temporal relationships in the MPB northern range data using the directed graph in Figure 12. Each set of nodes on the same value of the x-axis is a map of polygons, and edges are determined by the spatial relationships specified by the parameters of the stamp function, in this case with the distance threshold set to 2 km. The directed graph can be used to either analyze the structural properties of the evolution of the outbreak, or to identify specific events of interest. Here we queried the graph by degree to highlight nodes (i.e., polygons) with a high degree of interaction with polygons the next year, based on

Discussion
STAMP analysis enables a finer understanding of the spatial processes associated with moving polygons. Specifically STAMP facilitates the categorization of 11 different types of movement events (see RNBW), which provide information for evaluating the processes generating observed polygon movements. The STAMP categorizations are enhanced through distance and directional analysis, which adds geographical context to movement event analysis. For example, in the hurricane Katrina analysis we are able to clearly identify times and directions of movement when the hurricane was speeding up (high generation area in Figure 7) and slowing down (high disappearance area in Figure 7). These inferences are corroborated when compared with the Hausdorff distance metric (Figure 8), and would not have been observed using a traditional point-based analysis of hurricane trajectories. In the MPB example, we were able to compare directional change distributions, and isolate a specific event where infestation fragmentation occurred. These examples provide just a flavor of the types of analysis that can be accomplished with the STAMP framework as implemented in the stampr package.
There exist opportunities for new methods and techniques for polygon movement analysis. One of the main limitations of polygon movement analysis is that there is often no singular result with respect to either distance or directional relationships. The lack of a formal truth makes interpreting the results application specific, highlights the need for an exploratory data analysis approach, and requires analysts to carefully evaluate which methods align best with their study objectives. Future developments would benefit from considering null distributions of movement indices for hypothesized change types, which could then be used to assess the significance of observed events.
Within the stampr package we have included a range of methods for performing distance and direction analysis on polygon movements. The tools we have implemented represent the state-of-the-art for polygon movement analysis. However, we see opportunities for continued development of spatial-temporal polygon analysis methods. As new methods are documented we can add these to the package function set. Finally, we hope that the stampr package continues to build on the suite of spatial methods available within R, attracting a wide range users to polygon-based spatial analysis techniques.

Conclusions
Movement research in geographic information science is proliferating due to new datasets and the growing complexities of data (Long and Nelson 2013). The multi-method approach to polygon change analysis employed here is valuable in that it produces a broader understanding of the underlying movement dynamics, and facilitates comparisons between methods when dissimilar results are achieved. Here we highlight an example using the movement of hurricane Katrina and demonstrate the value of a polygon-based approach when contrasted with the more common trajectory-based analyses, and an analysis of landscape scale forest infestation spread. Future research should look to expand the application of polygon movement analysis now that these computational tools are readily available, as well as extending the methodology to incorporate spatial simulation models and other forms of spatial representation into the event-based framework implemented here. There are many application areas where polygon change analysis has been ignored (due to limited methods and tools) or where more simplistic approaches have been favored (based purely on polygon centroids). The development of a spatial toolkit for polygon change analysis, the stampr package, will reinvigorate research into performing spatial-temporal analysis of moving polygons.