# Important: If you haven't already, you need to install pygbif using the command below
# Run this command in the terminal to install: pip install pygbif

# Import packages
import os # For file paths, data management, interacting with the operating system
import pathlib # Object oriented file path management
import time # Time how long operations take in Python
import zipfile # Deal with zip files
from getpass import getpass # Security login for downloading GBIF data
from glob import glob # Search for specific file extensions

# Libraries for data analysis 
import pandas as pd # Data analysis for tabular data
import geopandas as gpd # Data analysis for geospatial data
import pygbif.occurrences as occ # Get occurrence data
import pygbif.species as species # Get species data
import requests # Download data

# Libraries for Dynamic mapping
import cartopy.crs as ccrs #Coordinate reference systems 
import panel as pn #customize the dynamic map
import hvplot.pandas #make the holoview plot with geo-dataframe
import calendar # Change month numbers to month names

# Define data directory in user's home directory
data_dir_2009 = os.path.join(
    # Home directory
    pathlib.Path.home(),
    # Earth analytics data directory
    'earth-analytics',
    'data',
    # Project directory
    'redstart-migration-2009',
)
print(data_dir_2009)
# Make the directory using the path defined above
os.makedirs(data_dir_2009, exist_ok=True)

C:\Users\warno\earth-analytics\data\redstart-migration-2009

# GBIF LOGIN CODE CHUNK
# Code credit: Elsa Culler
# This code allows us to download data in Python from gbif.org
# This code ASKS for your credentials for gbif.org
# and saves it for the rest of the session.
# NEVER put your credentials into your code

# GBIF needs a username, password, and email 
# All 3 need to match the account
reset = False

# Request and store username
if (not ('GBIF_USER'  in os.environ)) or reset:
    os.environ['GBIF_USER'] = input('GBIF username:')

# Securely request and store password
if (not ('GBIF_PWD'  in os.environ)) or reset:
    os.environ['GBIF_PWD'] = getpass('GBIF password:')
    
# Request and store account email address
if (not ('GBIF_EMAIL'  in os.environ)) or reset:
    os.environ['GBIF_EMAIL'] = input('GBIF email:')

# Search for the American Redstart (Setophaga ruticilla)
backbone = species.name_backbone(name='Setophaga ruticilla')

# Save the species key
species_key = backbone["speciesKey"]
backbone

{'usageKey': 2489985,
 'scientificName': 'Setophaga ruticilla (Linnaeus, 1758)',
 'canonicalName': 'Setophaga ruticilla',
 'rank': 'SPECIES',
 'status': 'ACCEPTED',
 'confidence': 99,
 'matchType': 'EXACT',
 'kingdom': 'Animalia',
 'phylum': 'Chordata',
 'order': 'Passeriformes',
 'family': 'Parulidae',
 'genus': 'Setophaga',
 'species': 'Setophaga ruticilla',
 'kingdomKey': 1,
 'phylumKey': 44,
 'classKey': 212,
 'orderKey': 729,
 'familyKey': 5263,
 'genusKey': 2489984,
 'speciesKey': 2489985,
 'class': 'Aves'}

gbif_pattern_2009 = os.path.join(data_dir_2009, '*.csv')
if not glob(gbif_pattern_2009):
    # Only submit one request
    if not 'GBIF_DOWNLOAD_KEY' in os.environ:
        # Submit query to GBIF
        gbif_query = occ.download([
            f"speciesKey = { species_key }",
            "hasCoordinate = True",
            "year = 2009",
        ])
        os.environ['GBIF_DOWNLOAD_KEY'] = gbif_query[0]
        print("Submitted query to GBIF")
    
    # Wait for download to build
    download_key = os.environ['GBIF_DOWNLOAD_KEY']
    wait = occ.download_meta(download_key)['status']
    while not  wait=='SUCCEEDED':
        wait = occ.download_meta(download_key)['status']
        time.sleep(5)
    print("Building download")

    # Download GBIF data
    print("Downloading GBIF Data")
    download_info = occ.download_get(
        os.environ['GBIF_DOWNLOAD_KEY'],
        path=data_dir_2009)
    print("Finished Downloading GBIF Data")
    
    # Unzip GBIF data
    with zipfile.ZipFile(download_info['path']) as download_zip:
        download_zip.extractall(path=data_dir_2009)
    print("Unzipping GBIF data")

# Find extracted .csv file path and take 1st result
gbif_path_2009 = glob(gbif_pattern_2009)[0]

with open(gbif_path_2009) as f:
    for i in range(2):
        print(f.readline().strip())

gbifID	datasetKey	occurrenceID	kingdom	phylum	class	order	family	genus	species	infraspecificEpithet	taxonRank	scientificName	verbatimScientificName	verbatimScientificNameAuthorship	countryCode	locality	stateProvince	occurrenceStatus	individualCount	publishingOrgKey	decimalLatitude	decimalLongitude	coordinateUncertaintyInMeters	coordinatePrecision	elevation	elevationAccuracy	depth	depthAccuracy	eventDate	day	month	year	taxonKey	speciesKey	basisOfRecord	institutionCode	collectionCode	catalogNumber	recordNumber	identifiedBy	dateIdentified	license	rightsHolder	recordedBy	typeStatus	establishmentMeans	lastInterpreted	mediaType	issue
985926909	4fa7b334-ce0d-4e88-aaae-2e0c138d049e	URN:catalog:CLO:EBIRD:OBS134163583	Animalia	Chordata	Aves	Passeriformes	Parulidae	Setophaga	Setophaga ruticilla		SPECIES	Setophaga ruticilla (Linnaeus, 1758)	Setophaga ruticilla		US	Boyne City	Michigan	PRESENT	1	e2e717bf-551a-4917-bdc9-4fa0f342c530	45.21671	-85.01395							2009-05-18	18	5	2009	2489985	2489985	HUMAN_OBSERVATION	CLO	EBIRD	OBS134163583				CC_BY_4_0		obsr180729			2025-11-05T07:42:44.891Z		CONTINENT_DERIVED_FROM_COORDINATES;TAXON_CONCEPT_ID_NOT_FOUND

# Load the GBIF data, this time defining the delimiter (tabular) represents as '\t'
# We also set the index to gbifID, and select the columns of interest (lat, long, and month)
gbif_df = pd.read_csv(
    gbif_path_2009,
    delimiter='\t',
    index_col='gbifID',
    usecols=['gbifID','decimalLongitude','decimalLatitude','month']
)
gbif_df.head()

# Converting the dataframe to a geodataframe
# Using geopandas library
gbif_gdf = (
    gpd.GeoDataFrame(
        gbif_df, 
        geometry=gpd.points_from_xy(
            gbif_df.decimalLongitude, 
            gbif_df.decimalLatitude), 
        crs="EPSG:4326") # Manually set the crs to match the ecoregions data
    # Select the desired columns
    [['month','geometry']]
)
gbif_gdf

# Ecoregions data from GeographyRealm, RESOLVE, and OneEarth: 
# GeographyRealm: https://www.geographyrealm.com/terrestrial-ecoregions-gis-data/
# OneEarth: https://www.oneearth.org/announcing-the-release-of-ecoregion-snapshots/
eco_url = (
    "https://storage.googleapis.com/teow2016/Ecoregions2017.zip")

# Set up a path to save the data to your machine
eco_dir = os.path.join(data_dir_2009, 'resolve_ecoregions') 
# Make ecoregions directory
os.makedirs(eco_dir, exist_ok=True)

# Join ecoregions shapefile path
eco_path = os.path.join(eco_dir, 'ecoregions.shp')

# Only download once
if not os.path.exists(eco_path):
    my_gdf = gpd.read_file(eco_url)
    my_gdf.to_file(eco_path)

# Opening the download and renaming some of the column names
eco_gdf = (
    gpd.read_file(eco_path)
    [['OBJECTID','ECO_NAME','SHAPE_AREA','geometry']]
    .rename(columns={
        'OBJECTID': 'ecoregion_id',
        'ECO_NAME': 'name',
        'SHAPE_AREA': 'area'
    })
    .set_index('ecoregion_id')
)

# Plot the ecoregions quickly to check download
eco_gdf.plot()
eco_gdf.head()

gbif_ecoregion_gdf = (
    eco_gdf
    # Match the CRS of the GBIF data and the ecoregions
    .to_crs(gbif_gdf.crs)
    # Find ecoregion for each observation
    .sjoin(
        gbif_gdf,
        how='inner', 
        predicate='contains')
)
gbif_ecoregion_gdf

# Count the observations in each ecoregion each month
occurrence_df = (
    gbif_ecoregion_gdf
    # Reset index
    .reset_index()
    # For each ecoregion, for each month...
    .groupby(['ecoregion_id','month'])
    # ...count the number of occurrences
    # creating a new column 'occurrences' that counts all the gbif IDs in a group
    .agg(
        occurrences=('gbifID', 'count'),
        area=('area','first')) 
)

#Normalize by area
occurrence_df['density'] = (
    occurrence_df.occurrences
    / occurrence_df.area
)

# Get rid of rare observations (possible misidentification?)
occurrence_df = occurrence_df[occurrence_df.occurrences>1]

# Take the mean by ecoregion
mean_occurrences_by_ecoregion = (
    occurrence_df
    .groupby('ecoregion_id')
    .mean()
)
mean_occurrences_by_ecoregion

# Take the mean by month
mean_occurrences_by_month = (
    occurrence_df
    .groupby('month')
    .mean()
)
mean_occurrences_by_month

# Normalize by space and time for sampling effort
occurrence_df['norm_occurrences'] = (
    occurrence_df[['density']]
    / mean_occurrences_by_ecoregion[['density']]
    / mean_occurrences_by_month[['density']]
)
occurrence_df

occurrence_df.norm_occurrences.plot.hist()

<Axes: ylabel='Frequency'>

# Plotting a scatter with regular "occurrences" by month
occurrence_df.reset_index().plot.scatter(
    x='month', y='occurrences', c='ecoregion_id', logy=True
)

<Axes: xlabel='month', ylabel='occurrences'>

# Plotting a scatter with normalized occurrences by month
occurrence_df.reset_index().plot.scatter(
    x='month', y='norm_occurrences', c='ecoregion_id', logy=True
)

<Axes: xlabel='month', ylabel='norm_occurrences'>

# Simplify the geometry to speed up processing
eco_gdf.geometry = eco_gdf.simplify(.1, preserve_topology=False)

# Change the CRS to Mercator for mapping
eco_gdf = eco_gdf.to_crs(ccrs.Mercator())

# Check that the plot runs in a reasonable amount of time
eco_gdf.hvplot(geo=True, crs=ccrs.Mercator())

# Join the occurrences with the plotting GeoDataFrame
occurrence_gdf = eco_gdf.join(occurrence_df[['norm_occurrences']])

# Remove any empty geometry rows. These can cause errors if they persist in the data. 
occurrence_gdf = occurrence_gdf[~occurrence_gdf.geometry.is_empty]

# Get the plot bounds so they don't change with the slider
# xmin, ymin, xmax, ymax = gbif_gdf.to_crs(ccrs.Mercator()).total_bounds

# This code above will plot the entire range of the data (even other continents/unusual sightings)
# I am just interested in the South-North American data, so I set the bounds manually.

xmin = -35000000.0
xmax = -2000000.0
ymax = 11000000.0
ymin = -4000000.0

# List month names on the slider
month_widget = pn.widgets.DiscreteSlider(
    options={
        calendar.month_name[month_num]: month_num
        for month_num in range(1,13)
    }
)

# Plot occurrence by ecoregion and month
migration_plot = (
    occurrence_gdf
    .hvplot(
        c='norm_occurrences',
        groupby='month',
        # Use background tiles
        geo=True, crs=ccrs.Mercator(), tiles='CartoLight',
        title="American Redstart Migration, 2009",
        xlim=(xmin, xmax), ylim=(ymin, ymax),
        frame_height=220,
        frame_width=220,
        widgets={'month': month_widget},
        widget_location='bottom'
    )
)

# Save the plot
migration_plot.save('redstart-migration-2009.html', embed=True)
# Display the plot
migration_plot

WARNING:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='p11486', ...)

BokehModel(combine_events=True, render_bundle={'docs_json': {'2ed70620-84ac-41e7-acaf-e70c900895bd': {'version…

occurrence_gdf.to_file("redstart_occurrence_2009.geojson", driver="GeoJSON")

INFO:Created 576 records

	decimalLatitude	decimalLongitude	month
gbifID
985926909	45.21671	-85.01395	5.0
985724600	37.55375	-77.46031	5.0
985687751	42.01461	-83.20015	9.0
985470319	42.76390	-70.80230	6.0
985469734	42.76390	-70.80230	6.0

	month	geometry
gbifID
985926909	5.0	POINT (-85.01395 45.21671)
985724600	5.0	POINT (-77.46031 37.55375)
985687751	9.0	POINT (-83.20015 42.01461)
985470319	6.0	POINT (-70.8023 42.7639)
985469734	6.0	POINT (-70.8023 42.7639)
...	...	...
1147062541	5.0	POINT (-76.3344 42.32292)
1147043863	5.0	POINT (-76.288 42.431)
1147037854	6.0	POINT (-73.35102 42.82333)
1038340085	5.0	POINT (-79.86786 38.65549)
1019440669	4.0	POINT (-75.7591 39.66956)

	name	area	geometry
ecoregion_id
1.0	Adelie Land tundra	0.038948	MULTIPOLYGON (((158.7141 -69.60657, 158.71264 ...
2.0	Admiralty Islands lowland rain forests	0.170599	MULTIPOLYGON (((147.28819 -2.57589, 147.2715 -...
3.0	Aegean and Western Turkey sclerophyllous and m...	13.844952	MULTIPOLYGON (((26.88659 35.32161, 26.88297 35...
4.0	Afghan Mountains semi-desert	1.355536	MULTIPOLYGON (((65.48655 34.71401, 65.52872 34...
5.0	Ahklun and Kilbuck Upland Tundra	8.196573	MULTIPOLYGON (((-160.26404 58.64097, -160.2673...

	name	area	geometry	gbifID	month
ecoregion_id
13.0	Alberta-British Columbia foothills forests	17.133639	MULTIPOLYGON (((-119.53979 55.81661, -119.5443...	3269620424	9.0
13.0	Alberta-British Columbia foothills forests	17.133639	MULTIPOLYGON (((-119.53979 55.81661, -119.5443...	1205984752	6.0
13.0	Alberta-British Columbia foothills forests	17.133639	MULTIPOLYGON (((-119.53979 55.81661, -119.5443...	5425832666	5.0
13.0	Alberta-British Columbia foothills forests	17.133639	MULTIPOLYGON (((-119.53979 55.81661, -119.5443...	5531404701	5.0
13.0	Alberta-British Columbia foothills forests	17.133639	MULTIPOLYGON (((-119.53979 55.81661, -119.5443...	5706956904	5.0
...	...	...	...	...	...
839.0	Northern Rockies conifer forests	35.905513	POLYGON ((-119.99977 54.53117, -119.8914 54.45...	256767436	7.0
839.0	Northern Rockies conifer forests	35.905513	POLYGON ((-119.99977 54.53117, -119.8914 54.45...	256438480	7.0
839.0	Northern Rockies conifer forests	35.905513	POLYGON ((-119.99977 54.53117, -119.8914 54.45...	256317445	6.0
839.0	Northern Rockies conifer forests	35.905513	POLYGON ((-119.99977 54.53117, -119.8914 54.45...	599681238	6.0
839.0	Northern Rockies conifer forests	35.905513	POLYGON ((-119.99977 54.53117, -119.8914 54.45...	599682307	6.0

	occurrences	area	density
ecoregion_id
13.0	6.666667	17.133639	0.389098
17.0	118.714286	7.958751	14.916196
23.0	3.600000	3.346216	1.075842
33.0	182.428571	16.637804	10.964703
34.0	104.714286	18.674884	5.607226
...	...	...	...
821.0	9.000000	26.034401	0.345696
827.0	2.000000	9.664680	0.206939
833.0	12.428571	0.610793	20.348238
838.0	6.000000	4.286144	1.399860
839.0	44.200000	35.905513	1.231009

	occurrences	area	density
month
1.0	8.414634	3.946549	34.251673
2.0	8.714286	5.638937	5.759934
3.0	8.756757	4.923813	8.343091
4.0	32.644444	10.337087	7.115234
5.0	183.621212	19.332582	17.410522
6.0	93.618182	23.795967	6.333724
7.0	51.365854	26.555008	3.139044
8.0	57.519231	22.281602	6.707274
9.0	102.031250	17.284795	33.416582
10.0	31.453125	10.457288	67.876760
11.0	9.971429	6.995604	38.831634
12.0	8.075000	4.866798	24.591023

GBIF Data Download and Methods for 2009¶