Range edges

Categorizing historical and modern range edges¶

In [1]:

import ecospat.ecospat as ecospat_full
from ecospat.stand_alone_functions import (
    get_species_code_if_exists,
    merge_touching_groups,
    assign_polygon_clusters,
    classify_range_edges,
    update_polygon_categories,
    get_start_year_from_species,
    fetch_gbif_data_with_historic,
    convert_to_gdf,
    process_gbif_data_pipeline,
    calculate_density,
    summarize_polygons_with_points,
    process_species_historical_range,
    analyze_species_distribution,
)

import ecospat.ecospat as ecospat_full from ecospat.stand_alone_functions import ( get_species_code_if_exists, merge_touching_groups, assign_polygon_clusters, classify_range_edges, update_polygon_categories, get_start_year_from_species, fetch_gbif_data_with_historic, convert_to_gdf, process_gbif_data_pipeline, calculate_density, summarize_polygons_with_points, process_species_historical_range, analyze_species_distribution, )

Step-by-step historical¶

In [2]:

# Range maps of over 600 North American tree species were created by Elbert L. Little, Jr. from 1971-1977
# First we need to load in the historical Little data for a tree species to an ecospat map

historic_map = ecospat_full.Map()
species_name = "Populus angustifolia"
code = get_species_code_if_exists(species_name)
historic_map.load_historic_data(species_name, add_to_map=True)
historic_map

# Range maps of over 600 North American tree species were created by Elbert L. Little, Jr. from 1971-1977 # First we need to load in the historical Little data for a tree species to an ecospat map historic_map = ecospat_full.Map() species_name = "Populus angustifolia" code = get_species_code_if_exists(species_name) historic_map.load_historic_data(species_name, add_to_map=True) historic_map

Out[2]:

In [3]:

# Next we need to remove lakes and major bodies of water and merge touching polygons

range_no_lakes = historic_map.remove_lakes(historic_map.gdfs[code])

# We can update the buffer_distance parameter based what polygons we want to merge; 5000m is a good start

merged_polygons = merge_touching_groups(range_no_lakes, buffer_distance=5000)

merged_polygons.plot()

# Next we need to remove lakes and major bodies of water and merge touching polygons range_no_lakes = historic_map.remove_lakes(historic_map.gdfs[code]) # We can update the buffer_distance parameter based what polygons we want to merge; 5000m is a good start merged_polygons = merge_touching_groups(range_no_lakes, buffer_distance=5000) merged_polygons.plot()

Out[3]:

<Axes: >

In [4]:

# Finally, we can classify the range edges of the historical range

# Identifies large core polygons
clustered_polygons, largest_polygons = assign_polygon_clusters(merged_polygons)

# Classifies range edges based on latitudinal and longitudinal position to core polygons
classified_polygons = classify_range_edges(clustered_polygons, largest_polygons)

# Updates polygon categories for polygons on islands
updated_polygon = update_polygon_categories(largest_polygons, classified_polygons)

updated_polygon.plot(column="category", legend=True, figsize=(10, 12))

# Finally, we can classify the range edges of the historical range # Identifies large core polygons clustered_polygons, largest_polygons = assign_polygon_clusters(merged_polygons) # Classifies range edges based on latitudinal and longitudinal position to core polygons classified_polygons = classify_range_edges(clustered_polygons, largest_polygons) # Updates polygon categories for polygons on islands updated_polygon = update_polygon_categories(largest_polygons, classified_polygons) updated_polygon.plot(column="category", legend=True, figsize=(10, 12))

No overlapping polygons found — returning original classifications.

Out[4]:

<Axes: >

In [5]:

# We can also plot these polygons on an ecospat map
historical_map_poly = ecospat_full.Map()
historical_map_poly.add_range_polygons(updated_polygon)
historical_map_poly

# We can also plot these polygons on an ecospat map historical_map_poly = ecospat_full.Map() historical_map_poly.add_range_polygons(updated_polygon) historical_map_poly

Out[5]:

Step-by-step modern¶

In [6]:

# First we need to fetch modern and historic GBIF data. Historic GBIF data will be used to calculate population density change.

# Let's retrieve the year the associated little map was published for this species
start_year = get_start_year_from_species(species_name)
start_year = int(start_year)

# Now we will pull 1000 GBIF occurrences from 2025 backwards and from 1976 (start year) backwards
data = fetch_gbif_data_with_historic(
    species_name,
    limit=1000,
    start_year=start_year,
    end_year=2025,
    continent="north_america",
)
modern_data = data["modern"]
historic_data = data["historic"]

# Finally, we convert this raw GBIF data into a gdf
historic_gdf = convert_to_gdf(historic_data)
modern_gdf = convert_to_gdf(modern_data)

# As an example, we will view the first few rows of the modern GBIF gdf
modern_gdf.head()

# First we need to fetch modern and historic GBIF data. Historic GBIF data will be used to calculate population density change. # Let's retrieve the year the associated little map was published for this species start_year = get_start_year_from_species(species_name) start_year = int(start_year) # Now we will pull 1000 GBIF occurrences from 2025 backwards and from 1976 (start year) backwards data = fetch_gbif_data_with_historic( species_name, limit=1000, start_year=start_year, end_year=2025, continent="north_america", ) modern_data = data["modern"] historic_data = data["historic"] # Finally, we convert this raw GBIF data into a gdf historic_gdf = convert_to_gdf(historic_data) modern_gdf = convert_to_gdf(modern_data) # As an example, we will view the first few rows of the modern GBIF gdf modern_gdf.head()

Out[6]:

	species	decimalLatitude	decimalLongitude	year	eventDate	basisOfRecord	geometry
0	Populus angustifolia	40.880712	-111.840163	2025	2025-01-05	HUMAN_OBSERVATION	POINT (-111.84016 40.88071)
1	Populus angustifolia	39.695967	-104.920282	2025	2025-01-09	HUMAN_OBSERVATION	POINT (-104.92028 39.69597)
2	Populus angustifolia	40.775008	-111.467794	2025	2025-02-07	HUMAN_OBSERVATION	POINT (-111.46779 40.77501)
3	Populus angustifolia	40.607400	-105.103214	2025	2025-03-02	HUMAN_OBSERVATION	POINT (-105.10321 40.6074)
4	Populus angustifolia	37.261355	-113.441772	2025	2025-03-05	HUMAN_OBSERVATION	POINT (-113.44177 37.26136)

In [7]:

import matplotlib.pyplot as plt

# Now we will need to processes this raw GBIF data in order to classify range edges

classified_modern = process_gbif_data_pipeline(
    modern_gdf,
    species_name=species_name,
    is_modern=True,
    end_year=2025,
    continent="north_america",
)

ax = classified_modern.plot(column="category", legend=True, figsize=(10, 12))
ax.set_title("Modern GBIF Range Edges")

classified_historic = process_gbif_data_pipeline(
    historic_gdf, is_modern=False, end_year=2025, continent="north_america"
)

ax_historic = classified_historic.plot(column="category", legend=True, figsize=(10, 10))
ax_historic.set_title("Historic GBIF Range Edges")

import matplotlib.pyplot as plt # Now we will need to processes this raw GBIF data in order to classify range edges classified_modern = process_gbif_data_pipeline( modern_gdf, species_name=species_name, is_modern=True, end_year=2025, continent="north_america", ) ax = classified_modern.plot(column="category", legend=True, figsize=(10, 12)) ax.set_title("Modern GBIF Range Edges") classified_historic = process_gbif_data_pipeline( historic_gdf, is_modern=False, end_year=2025, continent="north_america" ) ax_historic = classified_historic.plot(column="category", legend=True, figsize=(10, 10)) ax_historic.set_title("Historic GBIF Range Edges")

Out[7]:

Text(0.5, 1.0, 'Historic GBIF Range Edges')

In [8]:

# We then need to calculate the density of points (or unique individuals per polygon)

classified_modern = calculate_density(classified_modern)
classified_historic = calculate_density(classified_historic)

summarized_modern = summarize_polygons_with_points(classified_modern)

summarized_modern.head()

# We then need to calculate the density of points (or unique individuals per polygon) classified_modern = calculate_density(classified_modern) classified_historic = calculate_density(classified_historic) summarized_modern = summarize_polygons_with_points(classified_modern) summarized_modern.head()

Out[8]:

	geometry_id	geometry	category	AREA	cluster	n_points
0	1504486e4d36dd788c81f40012d15c4b	POLYGON ((-111.10275 44.76172, -111.08115 45.6...	leading (0.99)	7138.194457	1	16
1	17142a26759936a72cb3eb9ab97b1747	POLYGON ((-108.00279 45.98829, -108.52747 45.7...	leading (0.99)	2112.726453	1	12
2	1bb86ee22630c1e677e7cbf9fdd6afed	POLYGON ((-112.45636 34.55472, -112.42746 34.5...	relict (0.01 latitude)	124.508601	1	11
3	2686f6d8b0d0f7f01397443204c95b9a	POLYGON ((-114.22881 38.92158, -114.11357 39.3...	trailing (0.1)	1538.418666	1	9
4	2984dec013ddcad1ab7ad2b17a6b76f0	POLYGON ((-111.73118 35.15945, -111.57835 35.2...	relict (0.01 latitude)	194.434139	1	7

In [9]:

# Finally, lets add these modern polygons to an ecospat map

modern_map_poly = ecospat_full.Map()
modern_map_poly.add_range_polygons(summarized_modern)
modern_map_poly

# Finally, lets add these modern polygons to an ecospat map modern_map_poly = ecospat_full.Map() modern_map_poly.add_range_polygons(summarized_modern) modern_map_poly

Out[9]:

Pipeline functions¶

Historical pipeline¶

In [10]:

# Here we are going to generate the historic range map data
hist_pipeline = ecospat_full.Map()
hist_range = process_species_historical_range(
    new_map=hist_pipeline, species_name="Populus angustifolia"
)
hist_pipeline.add_range_polygons(hist_range)
hist_pipeline

# Here we are going to generate the historic range map data hist_pipeline = ecospat_full.Map() hist_range = process_species_historical_range( new_map=hist_pipeline, species_name="Populus angustifolia" ) hist_pipeline.add_range_polygons(hist_range) hist_pipeline

No overlapping polygons found — returning original classifications.

Out[10]:

Modern GBIF pipeline¶

In [11]:

classified_modern, classified_historic = analyze_species_distribution(
    "Populus angustifolia", record_limit=1000
)
classified_modern

classified_modern, classified_historic = analyze_species_distribution( "Populus angustifolia", record_limit=1000 ) classified_modern

Modern records (>= 1976): 1000
Historic records (< 1976): 254

Out[11]:

	point_geometry	year	eventDate	geometry	geometry_id	cluster	AREA	category	density
0	POINT (-111.840163 40.880712)	2025	2025-01-05	POLYGON ((-112.16175 40.79402, -112.1358 40.81...	cb869cb2640256dc655c5ffd650a009f	1	52264.281507	core	0.003597
1	POINT (-111.467794 40.775008)	2025	2025-02-07	POLYGON ((-112.16175 40.79402, -112.1358 40.81...	cb869cb2640256dc655c5ffd650a009f	1	52264.281507	core	0.003597
2	POINT (-110.869423 39.731437)	2025	2025-03-29	POLYGON ((-112.16175 40.79402, -112.1358 40.81...	cb869cb2640256dc655c5ffd650a009f	1	52264.281507	core	0.003597
3	POINT (-111.826781 40.765911)	2025	2025-04-27	POLYGON ((-112.16175 40.79402, -112.1358 40.81...	cb869cb2640256dc655c5ffd650a009f	1	52264.281507	core	0.003597
4	POINT (-111.830586 40.497927)	2025	2025-05-25	POLYGON ((-112.16175 40.79402, -112.1358 40.81...	cb869cb2640256dc655c5ffd650a009f	1	52264.281507	core	0.003597
...	...	...	...	...	...	...	...	...	...
724	POINT (-108.596663 45.713265)	2020	2020-10-04	POLYGON ((-108.00279 45.98829, -108.52747 45.7...	17142a26759936a72cb3eb9ab97b1747	1	2112.726453	leading (0.99)	0.005680
725	POINT (-108.59076 45.709948)	2019	2019-02-17	POLYGON ((-108.00279 45.98829, -108.52747 45.7...	17142a26759936a72cb3eb9ab97b1747	1	2112.726453	leading (0.99)	0.005680
726	POINT (-104.491531 37.019506)	2023	2023-06-03	POLYGON ((-104.46276 37.09835, -104.43328 36.8...	87e87496522b18430b3b3dcca4d8e97e	0	65.282891	trailing (0.05)	0.045954
727	POINT (-104.43328 36.885445)	2022	2022-06-13	POLYGON ((-104.46276 37.09835, -104.43328 36.8...	87e87496522b18430b3b3dcca4d8e97e	0	65.282891	trailing (0.05)	0.045954
728	POINT (-104.462763 37.098354)	2021	2021-09-12	POLYGON ((-104.46276 37.09835, -104.43328 36.8...	87e87496522b18430b3b3dcca4d8e97e	0	65.282891	trailing (0.05)	0.045954

729 rows × 9 columns

In [12]:

modern_pipeline_summary = summarize_polygons_with_points(classified_modern)
modern_pipeline_map = ecospat_full.Map()
modern_pipeline_map.add_range_polygons(modern_pipeline_summary)
modern_pipeline_map

modern_pipeline_summary = summarize_polygons_with_points(classified_modern) modern_pipeline_map = ecospat_full.Map() modern_pipeline_map.add_range_polygons(modern_pipeline_summary) modern_pipeline_map

Out[12]: