Point Pattern Analysis
For a theoretical discussion on point patter analysis, see the lecture notes.
Before tackling this tutorial, you will need to download and install a dataset following these instructions:
Create a folder called PPA somewhere under your personal directory (e.g. C:\Users\jdoe\Documents\Tutorials\PPA\).
Download the data for this exercise then extract the contents of PPA.zip into your newly created PPA folder.
Open the project file
Open the ppa.aprx file.
The map consists of two layers: A Beilschmiedia tree point vector layer (bei.shp) that shows the location of 3605 trees (it is assumed to be a complete census of trees); An elevation raster layer (dem.tif) showing elevation values above mean sea level (recorded in meters).
Uniform quadrat count
We’ll first create a uniform grid polygon vector layer.
Open the Create Fishnet tool from the Geoprocessing toolbox (Data Management >> Sampling >> Create Fishnet ).
Set the output file name to grid.shp.
Set the Template Extent to that of the dem.tif raster. This defines the study extent within which the trees were recorded.
Selecting the dem.tif raster in the last step will automatically populate the fishnet extent fields.
You can define the grid intervals by grid cell size or grid dimensions. Here, we'll opt for the latter approach.
Define the grid as being 5 rows tall and 10 columns wide.
Uncheck the Create Label Points box (Only keep this checked if you want the tool to generate a point layer representing the cell centroids).
Set the output type to Polygon.
Click Run.
The new grid shapefile should now be present in your Contents pane.
Next, we will tally the number of trees within each grid cell.
Select the grid layer in the Contents pane if it is not already selected.
In the Feature Layer tab, click on Data >> Add Spatial Join.
Make sure that grid is the target feature set.
Set bei point layer as the join feature.
Set the Match Option to Contains.
Other parameters can be left to their default settings.
Click OK.
The geoprocess added a Join_Count column that provides the number of points in each grid cell.
Because the cell sizes are the same, the modifiable areal unit problem is not an issue insofar that area is not influencing our perception of the point distribution. But changing cell sizes, even when kept uniform, can change the overall "look" of the density distribution.
Tessellated quadrat count
In the last step, you created quadrats from a uniform grid. Quadrats need not be uniform in shape and in size. In this step, you will create quadrats based on elevation intervals for the purpose of counting the number of trees at different elevation ranges. We'll first reclassify the continuous raster into a coded integer raster where each code represents an elevation range. We'll then convert the integer raster to a polygon vector layer that will be used to tally the points,
In the geoprocessing pane, search for Reclassify. Click on the second link (i.e. the Reclassify tool available in the Spatial Analyst toolbox).
In the Reclassify tool, select dem.tif as the input raster.
Make sure that VALUE is the selected reclassification field.
Under the reclassification table, click on the Classify link.
Set the classification method to Equal Interval.
Set the number of classes to 4.
Click OK.
The classification table should be automatically populated with values. You will keep the default reclass values of 1,2,3 and 4.
Name the output raster tess.tif (make sure not to save it in a geodatabase).
Click Run.
Next, we’ll convert the integer raster to a vector polygon layer. Note that this new raster has an attribute table whereas the original raster did not. This is to be expected since the newly created raster has just four unique integer values (ArcGIS will not create a raster attribute table if there are too many unique pixel values).
Note that you should not confuse the Count field in this new raster with the one we'll be creating from the point data. This count field tallies the number of pixels belonging to each class.
Next, we'll convert this raster to a polygon layer.
In the geoprocessing pane, search for the Raster to Polygon tool.
In the geoprocessing tool, set tess.tif as the input raster.
Make sure that the Value field is selected.
Set the output filename to quadrat_tess.shp.
Make sure to check the box Create multipart features.
Click Run.
By checking the multipart feature option, you ensure that each unique value appears just once in the attribute table. So if multiple polygons share the same attribute value, they will all point to the same record in the attribute table creating a one-to-many relationship between record and features.
Next, we’ll tally the number of points in each tessellated surface using the spatial join workflow described earlier.
Select the quadrat_tess.shp layer in the Contents pane if it is not already selected.
In the Feature Layer tab, click on Data >> Add Spatial Join.
Make sure that quadrat_tess is the target feature set.
Set bei as the joining feature.
Set the Match Option to Contains.
Other parameters can be left to their default settings.
Click OK.
The Join_Count attribute of the output layer stores the number of points falling in each elevation quadrats.
We could compare the counts between all four elevation intervals; however, the area of each interval may not be the same (recall the modifiable aerial unit problem). So we should normalize the count to area. This will require that we create a new attribute field: Density. But first, we will export this layer to a new layer called tess_join.shp to make the join permanent (this may be necessary when adding new fields on certain versions of the software).
Right-click the quadrat_tess layer, select Data >> Export Features.
Name the new file tess_join.shp.
Open the tess_join attribute table.
In the tess_join attribute table, click on the Add Field button .
This will bring up the Fields tab.
Create two new fields: Area and density and make sure to set them to float. (See the Editing attribute tables tutorial for a refresher on adding new fields)
Click Save in the Fields ribbon.
Right-click the Area field and select Calculate Geometry.
Select Area as the geometric property and Square Kilometers as the unit. Note that if the projection does not preserve area, you should select Area (geodesic) for the geometric property.
Click OK.
Right-click the Density field and select Calculate field.
Set the expression to !Join_Count! / !Area! .
Click OK.
Symbolize the layer by its density value.
Finally, we can make use of ArcGIS’ chart tool to plot count vs elevation range.
With the tess_join layer selected, open the Scatter Plot tool from the Data ribbon.
In the Chart Properties pane, map the gridcode variable to the x-axis and the Density variable to the y-axis.
Uncheck the Show linear trend options box.
Keep the other defaults.
You'll note that the color scheme adopted by the plot matches the color scheme used in the map. This facilitates comparison. Also, the figures are dynamically linked in that selecting an element in one figure will highlight the matching element(s) in the other.
The plot suggests a peak density at the third elevation interval. Referring back to the reclassification table, this matches the elevation range of 139.6 m to 149.6 m.
ArcGIS' point density tool
In this step we'll create a density map of tree counts using ArcGIS' Point Density tool.
In the Geoprocessing pane, open the Point Density tool from the Spatial Analyst Tools >> Density toolset. You can access the tool via the Toolbox tree (as shown in the accompanying figure), or by typing Point Density in the Search window.
Set the input feature field to bei.
Name the output raster point_density.tif.
Set the cell size to 10 (these are in mapping units which is meters in our example).
Define the neighborhood as a 3 pixel by 3 pixel rectangular window.
Define the output units in square kilometers. This will give us density values in counts per square kilometers.
Since the geoprocess will create a new raster from a vector layer, we need to explicitly define the raster extent since the tool will default to the smallest rectangle encompassing the input point layer.
Click on the Environments tab.
Set the Extent to that of the dem.tif layer.
Click Run.
Next, we'll modify the new raster layer's symbology.
Bring up the point_density.tif layer's Symbology.
In the Symbology pane, set the method to Equal Interval and the number of classes to 10.
Set the color scheme to a green hue.
We'll modify the first class to encompass just the 0 values. This will allow us to assign no color to zero valued pixels.
Change the first class' upper value to 0.
Click on the color swatch next to the upper value of 0 and set the color to no color.
Your output should look something like this.
Each pixel is assigned the number of points within a 3x3 pixel search window and that value is then divided by the area of that search window. For example, if a cell has one point inside the search window, its output value will be:
1 point / (9 cells * 0.0001 km2 area per cell) = 1111.1 points per square kilometer (within a 3x3 search window)
ArcGIS' kernel density tool
ArcMap offers two density tools: a point density tool (covered in the last section) and a kernel density tool. This would seem to suggest that these two tools perform distinct operations when in fact the point density tool is simply a special case of a kernel density function whereby all input values in the focal operation are assigned equal weight. ArcMap’s Kernel density tool applies a quartic function that assigns a non-uniform weight to each point based on its proximity to the output cell's center. This tool tends to generate smoother looking density rasters.
In the Geoprocessing pane, open the Kernel Density tool from the Spatial Analyst Tools >> Density toolset. You can access the tool via the Toolbox tree (as shown in the accompanying figure), or by typing Kernel Density in the Search window.
The Kernel density parameters differ slightly from those found in the Point Density tool. Most notably, it does not offer kernel options other than the default quartic kernel. It also limits the search window shape to a circle.
Another difference between this tool and the Point Density tool is the added geodesic distance option which you might find useful when analyzing large surface extents.
Point the input feature field to bei.
Name the output raster kernel_density.tif.
Set the cell size to 10 (these are in mapping units which is meters in our example).
Set the search radius to 15 map units.
Define the output units in square kilometers.
Keep the other parameter defaults. We'll stick with the Planar method given that this projection does a good job in preserving area/distance measurements. Note that if the coordinate system used does a poor job in preserving diatance, you should opt for the Geodesic Method.
You will also need to define the output extent to match that of the DEM layer under the Environments tab as was done with the Point density tool.
Click Run.
Next, generate a symbology scheme similar to the one used with the point density raster (make sure to set 0 pixel values to no color).
Nearest neighbor tool
So far, we’ve explored density based approaches to quantifying point patterns. Density based analysis usually focuses on a point pattern’s first order property—i.e. its distribution vis-à-vis location. Another property of interest is a point pattern’s spatial interaction, a second order effect. A statistic that can be used to quantify a point pattern’s second order property is the average nearest neighbor (ANN) statistic.
In your Geoprocessing pane, search for the Average Nearest Neighbor tool. You can also access the tool via Spatial Statistics Tools >> Analyzing Patterns .
Set bei as the input feature.
Keep the default Euclidean distance option.
Check the box next to Generate Report.
Set the study extent area to 500,000. The units are adopted from the layer's coordinate system--meters in this example.
Click Run.
The output is not a data layer but a report saved as an HTML file. The filename is listed in the Geoprocess' Details window.
In the Analysis tab, click on the History button.
Click on the Average Nearest Neighbor tool. This will bring up a side window. (You might need to hover your cursor above the tool for a second or two to see the window pop up).
In the side window, click on Report File link (this should open an HTML file saved in your project folder).
ArcGIS should open the file in a web browser.
The infographic at the top of the report suggests that the trees are more clustered than expected had the trees been distributed at random across the study extent. The p-value of <0.001 tells us that there is less than a 0.1 percent chance that we would be wrong in dismissing the null hypotheses (that being that the underlying process that generated the observed pattern is completely random).
The bottom of the report gives us the observed mean distance between trees in our study area. The ANN value is 4.33 meters. In other words, on average, the trees are 4.33 meters apart from their first nearest neighbor.
A word of caution: This tool does not take the shape of the study area into consideration. Hence, you should interpret the hypothesis test results with caution if the study extent's shape differs from a square or rectangle. See the lecture notes for more detail.