Creation of a Reusable Habitat Suitability Model for The Eastern Coyote:
Hamlet of Water Mill, Long Island New York

Abstract

The Eastern coyote (Canis latrans var.) has become ubiquitous throughout New York State, recently being documented breeding on Long Island, the last region of the state thought to be without them. Due to the behavioral ecology of the Eastern coyote and Long Island's lack of any relatively large predators, the region is primed for negative human-widlife interactions if proactive management does not begin; additionally, more empirical research is needed on the Eastern coyote in the region to better inform and assist. Using ArcGIS Pro, Python and procedures established in U.S Fish & Wildlife ESM 103 I created a predictive Habitat Suitability Model (HSM) to assist with future research and management efforts. The resulting model is easily shareable and reusable. The intended application is to assist future research efforts that seek to better understand coyote abundance and range, in addition to allowing the identification of priority areas for management. Here I showcase use of the model to produce a predictive surface in a pilot analysis within the Hamlet of Water Mill, Long Island NY. Watermill is a small hamlet on the south-east fork of Long Island with a recent confirmed coyote sighting. The hamlet also has a diverse landscape that is somewhat of a microcosm for the rest of the region, making it an ideal area of interest (AOI) to test the model.

Background

The Eastern coyote is larger than it's western cousin typically reaching weights of 30 to 40 pounds and heights up to 24 inches (Gompper, 2002). It was long suspected and more recentely confirmed that Eastern coyotes are the result of hybridization with wolves & dogs. Additionally, larger prey that is more readily available is another non-exlusive factor to the Eastern coyotes increased body size (Gompper 2002; Way & Lynn 2016). A common sight throughout much of New York, though up until the early 20th century (1930s - 1940s) they were mostly limited to the mid and south-western areas of the United States (Bogan, 2014; Nagy et al., 2017). Despite their numbers being estimated as 20,000 – 30,000 in the state they are relatively uncommon in the southern urban counties and Long Island (ESF 2006). Some breeding populations have been confirmed within Bronx county parks; but no breeding populations were thought to be on Long Island, until 2016 when a breeding pair with four pups was photographed with a camera trap in northern Queens county (Nagy et al., 2017). Unfortunately, the breeding pair and most of the pups were euthanized by authorities despite public opinion largely supporting coexistence. Despite this, researchers are confident that these coyotes had likely been breeding for several years, there are more breeding pairs on the island, and their expansion is inevitable (Nagy et al., 2017). The expansion of the coyote could potentially be beneficial for the island, where local ecosystems lack top down control from a predator, which can lead to the overabundance of some species (e.g. deer) which in turn negatively impact the ecosystems in which they inhabit. For example, deer over grazing reduces resource availability for other species and can lead to their decline (Bush et al., 2012; Shelton et al., 2014). This was seen with ravens in northern New York, which were all but removed from the area after the departure of the wolf. With the establishment of the coyote to the area, raven populations have recovered (ESF 2006). Inversely however, it has been shown that coyotes can limit of the abundance of other mid to small size predators, such as the red fox, which has important implications in terms of their effect on incidences of Lyme disease (Levi et al., 2012); though this is still a matter of debate (Way and White, 2013). Regardless, due to the almost certainty that we will see an expansion of the coyote population on Long Island coupled with the fact they are notoriously difficult to cull (Nagy et al., 2017), it would be best to develop proactive management strategies to maximize benefits, minimize negative impacts, and promote coexistence.

To this end, a HSM can assist in identifying priority areas which has numerous benefits. For example, citizens can know where they're are likely to encounter coyotes and prepare, and we will also know where educational efforts and signage for coyotes needs to be focused (Bogan, 2014; Nagy et al., 2017). Another benefit of a HSM is as a starting point for future research, such as abundance studies that will need to decide on the best areas to lay out transects for collecting samples and/or setting camera traps. The Hamlet of Water Mill made the best choice for an area of interst on the island due to it's' recent Eastern coyote sightings and diverse landscape. The hamlet is a good mix of forested areas to the north, developed in the south/south east, wetlands dispersed throughout with agricultural land in the center; this mix of land-cover is typical throughout much of the northeastern United States.

What is a Habitat Suitability Model?

A type of suitability model that assesses the quality of habitat for a species of interest. An HSM consists of components that represent a need of the species, they are broken down further into habitat variables, environmental variables that if modifeid are thought to affect the ability of the habitat to support the species of interest. The model results in a Habitat Suitability Index (HSI), an arbitrary scoring within an evaluation scale set by the maker of the model (U.S Fish & Wildlife Service ESM 103, 1981).

I identified 3 Main components for this model through review of the relevant literature. Each component is described via habitat variables:

  1. Resource Availability
    • Forest Edge
    • Canopy Cover
  2. Cover & Reproduction
    • Land Cover
  3. Anthropogenic Disturbance
    • Distance to Roads
    • Building Density

Figure 1: Conceptual Model of the HSM

Conceptual Model of the HSM

Methodology: Utilizing ArcPy to Create the Model

The structure of the model is such that each factor is calculated separately and represented by a raster layer, with values from some arbitrary evaluation scale (1 - 100 was chosen here). The lower a score, the lower the expected carrying capacity of the area and therefore less suitable the habitat. After each factor is calculated the rasters are added together and then divided by the number of total number of rasters (i.e. arithmetic mean). There are multiple ways to calculate a HSI, arithmetic mean is the simplest and was determined to be best suited for this situation due to the strong compensatory relationship that likely exists between the different habitat variables due to the generalist and oppurtunistic nature of the Eastern coyote (U.S Fish & Wildlife Service ESM 103, 1981).

To start the needed modules are imported and the working environment is set up, we need both ArcPy and the ArcPy spatial analysis module (for working with rasters).

In [1]:
import arcpy
from arcpy.sa import *

# set up environment
wks = r"coy_hsi_calculator\sample_data\analysis_output.gdb"
aoi = r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_boundary"
arcpy.env.overwriteOutput = True
arcpy.env.extent = (aoi)

Before any modeling can begin, a land-cover raster must be sourced with acceptable accuracy and spatial resolution for the intended purpose. For this study a custom land-cover raster was created with NAIP 1M Imagery using the NLCD 2011 land-cover classifcation system (81% Overall accuracy, 76% Kappa). The creation of the land-cover raster is beyond the scope of this analysis, but a video explanation can be found here. Below is a map showcasing the land cover for the Hamlet of Water Mill

Figure 2: Land Cover Map of the Hamlet of Water Mill

Land Cover Map

Land Cover Factor

The simplest factor to create is the land-cover factor. This factor represents cover and reproduction within the model. Eastern coyotes have been shown to have home ranges centering around forested areas during breeding, gestration and pup-rearing seasons with most den sites noted to be in open hardwoods or in small clearing surrounded by hardwoods (Pearson and Hirth 1991). In terms of denning and breeding most evidence indicates that disturbed, deciduous forest areas with open canopies offer the best cover for the Eastern coyote (Pearson and Hirth, 1991; Way et al., 2004; Kays et al., 2008; Nagy et al., 2017).

To create this factor, first I bring it into the workspace and then it is reclassified to assign scores from the evaluation scale to each class. A simple wrapper function is created that wraps of the ArcPy reclassify tool.

In [3]:
lc_ras = arcpy.Raster(r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_lc_2017")
lc_remap = r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_lc_remap_table"

# remap values:
# 1  (Shadow) - 1
# 11 (Water) - 1
# 20 (Developed) - 10
# 21 (Developed, Open Space) - 40
# 30 (Barren Land) - 25
# 41 (Deciduous Forest) - 100
# 42 (Evergreen Forest) - 85
# 80 (Agricultural/Crop/Pasture) - 70
# 90 (Wetlands) - 55

def reclassify_wrapper(raster,remap_obj,field = "Value"):
    """
    Wrapper around the reclassify tool
    
    Args:
        raster: The raster to be reclassified
        remap_obj: a RemapValue or RemapRange object
        field: Field to act on, "Value" by default.
        
    Returns:
        A reclassified raster layer
  
    """

    reclass = Reclassify(raster,field,remap_obj)

    return reclass

lc_factor = reclassify_wrapper(lc_ras,lc_remap)

Canopy Factor

The next factor to calculate is canopy cover, one of the habtiat variables relating to resource availability. Thicker canopies tend to have less understory vegetative cover and less small mammals that coyotes typically hunt such as hares, voles, etc. Additionally, a thick understory is also likely to have more wild berries and other vegatative resources that coyotes take advantage of more frequently during the summer (Oostring 1956, Severinghaus and Brown 1956, Bittner and Rongstad 1982, Litvaitis et al., 1985 - cited by Kays et al., 2008; Crimmins et al., 2012). Overall, coyotes have been shown to be more abundant in open canopy areas with less woody structure and many natural edges (Kays et al., 2008). Canopy data was sourced from MRLC and has a spatial resolution of 30m. 30m is too coarse for this analysis but being the only available canopy cover dataset it was resampled to a resolution of 1m. Once resampled the next step is to remap the raster values to be on the chosen evaluation scale. The workflow is identical to the process for creating the land-cover factor

In [4]:
# Reclassify the canopy raster
canopy_ras = arcpy.Raster(r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_canopy")
canopy_remap = r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_canopy_remap_table"
# remap values:
# 0% to 10% - 50
# 11% to 30% - 100
# 31% to 50% - 75
# 51% to 70% - 50
# 71% to 100% - 25
canopy_factor = reclassify_wrapper(canopy_ras,canopy_remap)
In [10]:
# View Canopy Factor Result (darker areas have lower percent canopy cover)
canopy_factor
Out[10]:

Forest Edge Factors

In addition to canopy cover, forest edge was also identified as a habitat variable predictive of resource availability. Here I focus on two such environments, forested areas that edge agricultural areas, and forested areas that edge wetlands. Forest that edge agricultural areas was included due to the high resource availability for coyotes in agricultural areas coupled with higher mortality (Crete et al., 2001; Richer et al., 2002; Kays et al., 2008; Mitchell et al., 2015; Hinton et al., 2015). An interesting dichotomy, with coyotes having a choice of higher risk but greater reward. Coyotes have also been shown to be adept at taking advantage of small habitat patches within otherwise dangerous areas such as urban areas, a similar high risk/high reward situation (Gehrt et al., 2009). In contrast to agricultural areas, forested areas provide better cover and breeding habitat than agricultural areas, typically with with adequate resource availability and lower mortality (depending on the forest type and structure). Considering that coyote abundance is mostly driven by resource availability, I theorize that coyotes will likely have higher abundance in forested areas that closely edge an agricultural area. The idea being that they can range into the agricultural areas when it is safe and retreat quickly back into forested areas when it is not (Crete et al., 2001; Gompper 2002; Kays et al.,2008). The Eastern coyote will also likely take advantage of any human resources found, with some some evidence to show that coyotes in this type of environment tend to have a higher intake of "unknown" food sources (likely refuse) (Crete et al., 2001). Forested areas that edge wetlands offer different oppurtunities. For example, wetland and other natural edges may assist the Eastern coyote in hunting larger prey such as White-tailed deer. Deer are potentially more vulnerable near frozen wetland areas and forest edges in the winter, when they become more important to the Eastern coyotes diet (Patterson and Messier 2003; Crimmins et al., 2002; Kays et al., 2008). Kays et al., 2008 suggest that disturbed forested areas with natural edges, notably wetlands, likely serve as better cover from human hunting presence than more open areas.

A function was created that calcules any edge area; this function is generalized and could potentially be used outside of this analysis. After calculating the edge areas the outputs are reclassified to the evaluation scale.

In [6]:
def calc_edge_area(outer_raster, inner_raster, radius=200):
    """
    Calculates an "edge" raster; an area of interest (e.g. forest) 
    that edges a cover type of interest (e.g. wetlands, streets, etc). 
    
    Args:
        outer_raster: The habitat of interest that will make up the edge 
                      environment. This raster should be boolean, with
                      cells that are of the cover type of interest
                      (e.g. forest) having a score of 1 and all other cells
                      scored 0
        inner_raster: The inner habitat that the outer raster edges. Should
                      be a boolean raster where cells of interest are scored 1
                      and all others 0
        radius: Controls the width of the edge raster returned, by default 
                200 map units will be used.
                
    Returns:
        A raster of the edge area defined per the provided arguments
        
    """

    # Run focal stats tool to find cells that "edge" the input cover_type
    # output raster will look similar to the original, but cells that edge the
    # cover will now have a score of 1 instead of 0                                  
    inner_edge = FocalStatistics(inner_raster,
                                        NbrRectangle(3, 3, "CELL"),
                                        "MAXIMUM",
                                        True)

    # find cells in the outer edge raster that are adjacent the 
    # inner edge
    outer_cells_adj_inner = ((outer_raster == 1)&(inner_edge > 0))

    # Run focal stats to find cells within x map unit radius
    # of the boundary betweent he two habitat rasters
    nbr_circle = NbrCircle(radius,"MAP")
    within_x_of_edge = FocalStatistics(outer_cells_adj_inner,
                                       nbr_circle,
                                       statistics_type = "MAXIMUM",
                                       ignore_nodata = True)


    # Create new raster where cells of the outer habitat are of the cover type
    # of interest (e.g. forest) and within x radius of the habitat boundary
    edge_raster = ((outer_raster == 1)&(within_x_of_edge > 0))

    return edge_raster
In [7]:
# Bring in Forest, Wetland, Agricultural Rasters
forest_ras = arcpy.Raster(r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_forest_ras")
wetland_ras = arcpy.Raster(r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_wetland_ras")
ag_ras = arcpy.Raster(r"coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_ag_ras")
In [8]:
# View the forest raster
forest_ras
Out[8]:
In [9]:
# View the wetland raster
wetland_ras
Out[9]:
In [11]:
# Run calc_edge_area to find the edge boundary of forest and wetland, keep default of 200 map units (ft here)
forest_edge_wetland = calc_edge_area(forest_ras,wetland_ras)
In [12]:
# Viewing the result from the above, we can see all the forested within 200ft of wetlands
# The resulting raster gives these areas a score of 1 while any other area is a score of 0
forest_edge_wetland
Out[12]:
In [13]:
forest_edge_ag = calc_edge_area(forest_ras,ag_ras)
# reclassify the edge areas
# scores of 1 (edge area) scored 100, scores of 0 (non-edge) scored 1 
edge_remap = RemapValue([[1,100],[0,1]])
for_edge_wet_factor = reclassify_wrapper(forest_edge_wetland,edge_remap)
for_edge_ag_factor = reclassify_wrapper(forest_edge_ag,edge_remap)

Building and Street Factors

Lastly, factors are calculated for anthropogenic disturbances. It has been shown coyote abundance is highest further from roads in areas of low housing density (Kays et al., 2008). Interstingly, human disturbance has also been shown to effect coyote active hours, with them having increased nocturnal activity when humans are least active (Way et al., 2004). Despite their abundance being negatively assoicated with human activity, Eastern coyotes can still be found within close proximity to human land uses (Mitchell et al., 2015).

Two functions were developed to calculate these factors, each utilize the ArcPy rescaling tool as the literature could not provide enough evidence to define categories for scoring. All we can say is that as distance to roads/streets decreases (habitat is closer to them) and building density increases (amount of buildings per a user defined area) we expect coyotes to be less abundant. The Rescaling tool works by utilizing a mathematical function to assign suitability values to an input raster on a continuous scale (ESRI Docs). The user choses a function that best predicates the relationship between the phenomenon being modeled and the variable to which the phenomenon responds to. Here it was decided that the function that best represented the factors was the simple "small" and "large" functions. Each of these functions has a midpoint and spread. Scoring changes once the midpoint is reached, with the change being function dependant (see code docs below), the spread controls how quickly the values will increase or decrease as they move from the midpoint.

The Small function is used when smaller input values are more preferred. In contrast, Large is used when larger input values are more preferred.

In [14]:
buildings = r"C:\Users\cvric\OneDrive\Desktop\coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_buildings"
streets = r"C:\Users\cvric\OneDrive\Desktop\coy_hsi_calculator\sample_data\water_mill_data.gdb\wm_streets"

def rescale_wrapper(dist_layer,function,midpoint="#",
                    spread="#",scale_start=1,scale_end=100):
    """
    Wrapper around the rescale tool that limits the transformation
    functions used to "small" and "large". Only dist_layer and function
    are required; midpoint will default to the mean of the dist_layer, and
    spread to a value of 5.
    
    Args:
        dist_layer: The a raster layer that represents a disturbance factor.
        function: Transformation function to use (must be "small" or "large").
        midpoint: Values greater than midpoint decrease in preference("small").
                  Values greater than midpoint increase in preference("large").
        spread: How quickly values increase and decrease as they move from the
                midpoint.
        scale_start: The start of the evaluation scale (e.g. 1).
        scale_end: The end of the evaluation scale (e.g. 100).
    
    Returns:
        a rescaled raster layer
        
    Raises:
        Exception: Raised if non-valid function parameter is supplied
    
    """

    if function.lower() == "small":
        transformation_function = TfSmall(midpoint, spread)
    elif function.lower() == "large":
        transformation_function = TfLarge(midpoint, spread)
    else:
        raise Exception("Must use either the 'small' or 'large' function, "
                        "check documentation and try again")

    reclass = RescaleByFunction(dist_layer, transformation_function, scale_start, scale_end)

    return reclass

def calc_building_factor(buildings,cellsize,radius,midpoint="#",
                         spread="#",scale_start=1,scale_end=100):
    """
    Calculates the building factor for the Eastern Coyote Habitat Suitability
    Model. Building points per the specified radius are calculated to gauge
    building density and the a rescale is performed.
    
    Args:
        buildings: Point or polygon features that represent buildings.
        radius: Radius used to calculate building density per unit area
        midpoint: Values greater than midpoint decrease in preference.
        spread: How quickly values increase and decrease as they move from the
                midpoint.
        scale_start: The start of the evaluation scale (e.g. 1).
        scale_end: The end of the evaluation scale (e.g. 100).
    
    Returns:
        The building factor for the model per the specified arguments
        
    Raises:
        Exception: Raised if a feature class of the wrong shapetype is provided
    
    """
    shape_type = arcpy.Describe(buildings).shapeType
    if shape_type == 'Polygon':
        building_pts = arcpy.FeatureToPoint_management(buildings,
                                                       r"in_memory\buildPts",
                                                       "INSIDE")
    elif shape_type == 'Point':
        building_pts = buildings

    else:
        raise Exception("buildings must be polygon or point!")


    pt_density = PointDensity(building_pts, "NONE",cellsize,
                              NbrCircle(radius,"MAP"))

    factor = rescale_wrapper(pt_density,"small",midpoint,
                             spread,scale_start,scale_end)

    arcpy.Delete_management(r"in_memory\buildPts")

    return factor


def calc_street_factor(streets,cellsize,midpoint="#",
                       spread="#",scale_start=1,scale_end=100):
    """
    Calculates the street factor for the Eastern Coyote Habitat Suitability
    Model.Euclidean distance is calculated for cells in the raster, further
    a cell is from a street the higher it will be scored
    
    Args:
        streets: line feature representing streets/roads.
        cellsize: Cell size to use for euclidean distance calculation
        midpoint: Values greater than midpoint increase in preference.
        spread: How quickly values increase and decrease as they move from the
                midpoint.
        scale_start: The start of the evaluation scale (e.g. 1).
        scale_end: The end of the evaluation scale (e.g. 100).
    
    Returns:
        The building factor for the model per the specified arguments
        
    """

    dt_streets = EucDistance(streets,cell_size=cellsize)

    factor = rescale_wrapper(dt_streets,"large",midpoint,
                             spread,scale_start,scale_end)

    return factor

For the building factor, the default values are used for the parameters relative to the rescaling, which was found to deliver a reasonable result. A radius of 8200 map units (feet here) is used (yields an area of ~ 20km^2); this is a reasonable assumption for an Eastern coyote home range within the AOI. The land cover raster cell sized is used here (cell size should be continuious throughout the model), the defaults are kept for the evaluation scale parameters.

In [15]:
building_factor = calc_building_factor(buildings,lc_ras,8200)
In [16]:
# View building factor
# Darker areas are more building dense
building_factor
Out[16]:

The street factor is calculated using a midpoint of 200 and spread of 5. The closer a cell is to a street, the lower it's evaluation score.

In [17]:
street_factor = calc_street_factor(streets,lc_ras,200,5)
In [18]:
# View street factor
street_factor
Out[18]:

Putting It Altogether

The last step is to calculate the final HSI and then clip the result to the area of interest. Before the score is calculated the Con and IsNull tools are used to ensure that each factor raster edge remains intact, converting any no data values to 0. The result is a predictive surface clipped to the AOI.

In [19]:
factors = [for_edge_wet_factor,for_edge_ag_factor,building_factor,
               street_factor,lc_factor,canopy_factor]

# Ensures that resulting raster edges remain intact
# converts no data values to 0
factors = [Con(IsNull(factor),0,factor) for factor in factors]

hsi = sum(factors)/6

# clip the resulting hsi raster to the area of interest
desc = arcpy.Describe(aoi).extent
extent = "{} {} {} {}".format(desc.XMin,desc.YMin,desc.XMax,desc.YMax)

# The resulting model is written to the output workspace with the name "coyote_hsi"
coyote_hsm = arcpy.Clip_management(hsi,extent,wks+r"\coyote_hsi",aoi,"0","ClippingGeometry","MAINTAIN_EXTENT")

Conclusion: Viewing The Result And More!

Below is a map showcasing the resulting model. We can see that the highest scores are disturbed forested areas that edge agricultural and/or wetlands and are far from anthropogenic disturbances. More interesting than the result is the models reusability and the prospect of adjustments as more information becomes available. The hope is that an expert/researcher utilizes the model and makes adjustments as they see fit. For example, the model could be used to decide where the best location is for a local abundance study, the results of which could be used to calibrate the model closer to the ideal.

Figure 3: Map Showcasing Result of HSM

Map Showcasing Result of HSM

Python Script Tool and Contact Information

To help facilitate this further, python script tool has been created and is available to download freely here. I am also available for contact at cvricella2@outlook.com and welcome any questions regarding the model!

References

Bogan, D. (2014, June). Rise of the Eastern Coyote. New York State Conservationist.

Bush, E., Buesching, C., Slade, E. and Macdonald, D. (2012). Woodland Recovery after Suppression of Deer: Cascade effects for Small Mammals, Wood Mice (Apodemus sylvaticus) and Bank Voles (Myodes glareolus). PLoS ONE, 7(2), p.e31404.

Communications, E. (2006). The Coyote in New York State. [online] Esf.edu. Available at: https://www.esf.edu/pubprog/brochure/coyote/coyote.htm [Accessed 14 May 2019].

Crête, M., Ouellet, J., Tremblay, J. and Arsenault, R. (2001). Suitability of the forest landscape for coyotes in northeastern North America and its implications for coexistence with other carnivores. Écoscience, 8(3), pp.311-319.

Gehrt, S., Anchor, C. and White, L. (2009). Home Range and Landscape Use of Coyotes in a Metropolitan Landscape: Conflict or Coexistence?. Journal of Mammalogy, 90(5), pp.1045- 1057.

Gompper, M. (2002). The Ecology of Northeast Coyotes: Current Knowledge and Priorities for Future Research. Wildlife Conservation Society, (Paper No. 17).

Habitat evaluation procedures (HEP). (1981). Washington, D.C.: Division of Ecological Services, U.S. Fish and Wildlife Service, Dept. of the Interior, p.ESM 103.

Hinton, J., van Manen, F. and Chamberlain, M. (2015). Space Use and Habitat Selection by Resident and Transient Coyotes (Canis latrans). PLOS ONE, 10(7), p.e0132203.

Kays, R., Gompper, M. and Ray, J. (2008). LANDSCAPE ECOLOGY OF EASTERN COYOTES BASED ON LARGE-SCALE ESTIMATES OF ABUNDANCE. Ecological Applications, 18(4), pp.1014-1027.

Levi, T., Kilpatrick, A., Mangel, M. and Wilmers, C. (2012). Deer, predators, and the emergence of Lyme disease. Proceedings of the National Academy of Sciences, 109(27), pp.10942-10947.

Mitchell, N., Strohbach, M. W., Pratt, R., Finn, W. C., & Strauss, E. G. (2015). Space use by resident and transient coyotes in an urban–rural landscape mosaic. Wildlife Research, 42(6), 461. doi:10.1071/wr15020

Nagy, C., Weckel, M., Monzón, J., Duncan, N. and Rosenthal, M. (2017). Initial colonization of Long Island, New York by the eastern coyote, Canis latrans (Carnivora, Canidae), including first record of breeding. Check List, 13(6), pp.901-907.

Person, D. K., & Hirth, D. H. (1991). Home Range and Habitat Use of Coyotes in a Farm Region of Vermont. The Journal of Wildlife Management, 55(3), 433. doi:10.2307/3808971

Richer, M.-C., Crête, M., Ouellet, J. P., Rivest, L. P., and Huot, J. (2002). Thelow performance of forest versus rural coyotes in northeastern NorthAmerica: inequality between presence and availability of prey. Ecoscience 9, 44–54.

Shelton, A., Henning, J., Schultz, P. and Clay, K. (2014). Effects of abundant white-tailed deer on vegetation, animals, mycorrhizal fungi, and soils. Forest Ecology and Management, 320, pp.39- 49.

Way, J. and White, B. (2013). Coyotes, Red Foxes, and the Prevalence of Lyme Disease. Northeastern Naturalist, 20(4), pp.655-665.