Density vs. Difference#
We want determine if the number of observations per pixel in the iNaturalis map relates to the difference between the sPlotOpen and iNaturalist map.
This section covers:
Load Data
Count observations per grid cell
Observation density vs. Discrepancy between iNat and sPlotOpen
# packages
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.ticker as ticker
Load Data#
Load iNaturalist-TRY data and sPlotOpen community weighted means
iNat_TRY = pd.read_csv("iNat_TRY_log.csv")
iNat_TRY.head(2)
| gbifID | scientificName | decimalLatitude | decimalLongitude | eventDate | dateIdentified | Dispersal unit length | Leaf Area | SLA | Leaf C | ... | Leaf delta15N | Leaf N P ratio | Leaf P | Plant Height | Seed mass | Seed length | Seeds per rep. unit | Stem conduit density | SSD | Conduit element length | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1229615436 | Commelina communis | 35.987483 | -79.057546 | 2013-07-07T00:00:00 | 2013-07-07T20:33:11 | NaN | NaN | NaN | NaN | ... | NaN | 2.5362 | 0.536493 | NaN | 2.13771 | NaN | NaN | NaN | NaN | NaN |
| 1 | 3384000233 | Commelina communis | 42.093762 | -75.923660 | 2021-08-23T13:06:06 | 2021-09-17T21:15:37 | NaN | NaN | NaN | NaN | ... | NaN | 2.5362 | 0.536493 | NaN | 2.13771 | NaN | NaN | NaN | NaN | NaN |
2 rows × 24 columns
sPlot = pd.read_csv("sPlotOpen/cwm_loc.csv")
sPlot.head(2)
| PlotObservationID | TraitCoverage_cover | Species_richness | TraitCoverage_pa | Leaf Area | SSD | SLA | Leaf C | Leaf N per mass | Leaf P | ... | Height_shrubs_lowest | Height_herbs_average | Height_herbs_lowest | Height_herbs_highest | SoilClim_PC1 | SoilClim_PC2 | Resample_1 | Resample_2 | Resample_3 | Resample_1_consensus | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 16 | 0.277778 | 3 | 0.333333 | 3.678311 | -1.047293 | 2.890748 | 6.128157 | 2.873263 | 1.114036 | ... | NaN | NaN | NaN | NaN | -3.66 | 0.546 | True | False | False | True |
| 1 | 17 | 0.038462 | 2 | 0.500000 | 3.678311 | -1.047293 | 2.890748 | 6.128157 | 2.873263 | 1.114036 | ... | NaN | NaN | NaN | NaN | -3.66 | 0.546 | True | False | False | True |
2 rows × 86 columns
Count observations per grid cell#
def global_grid_count(df, long, lat, deg, variables):
# create new dataframe to save the average value of each grid cell and variable
grouped_df = pd.DataFrame()
# convert degree into step size
step = int((360/deg) + 1)
bins_x = np.linspace(-180,180,step)
bins_y= np.linspace(-90,90,int(((step - 1)/2)+1))
# group latitude and longitude coordinates into bins
# create new columns 'x_bin' and 'y_bin'
df['x_bin'] = pd.cut(df[long], bins=bins_x)
df['y_bin'] = pd.cut(df[lat], bins=bins_y)
# raster coordinates are in center of raster cell
df['x_bin'] = df['x_bin'].apply(lambda x: ((x.left + x.right) /2) )
df['y_bin'] = df['y_bin'].apply(lambda x: ((x.left + x.right) /2) )
grouped_df = df.drop_duplicates(subset=['x_bin', 'y_bin'], keep='last')
grouped_df = grouped_df[['x_bin', 'y_bin']]
for v in variables:
sub_df = df[['y_bin', 'x_bin', v]]
# for each lat/lon group get count
grouped_v = sub_df.groupby(['x_bin', 'y_bin'], as_index=False)[v].count()
grouped_df = pd.merge(grouped_df, grouped_v,
on= ['x_bin', 'y_bin'],
how='left')
return grouped_df
# get counts
deg = [2]
trait = iNat_TRY.columns[6:24]
for d in deg:
df_iNat = global_grid_count(iNat_TRY, 'decimalLongitude', 'decimalLatitude', d, trait)
df_sPlot = global_grid_count(sPlot, 'Longitude', 'Latitude', d, trait)
# reshape data, so that we have only one Trait column
df_iNat_t = df_iNat.melt(id_vars=["x_bin", "y_bin"],
value_name="Count_iNat",
var_name="Trait")
df_sPlot_t = df_sPlot.melt(id_vars=["x_bin", "y_bin"],
value_name="Count_sPlot",
var_name="Trait")
# merge sPlot and iNat data into one dataframe
df_merged = pd.merge(df_sPlot_t, df_iNat_t, on=["x_bin", "y_bin", "Trait"] )
# keep only lines where we have a pixel in both datasets
df_merged = df_merged.dropna()
# save result to csv
filename="grid_count_" + str(d) + "_deg.csv"
df_merged.to_csv(filename, index=False)
Compare observation density per grid cell to scaled difference between sPlot and iNaturalist#
Open file that saved the counts per grid cell:
filename="grid_count_2_deg.csv"
counts = pd.read_csv(filename)
counts.head()
| x_bin | y_bin | Trait | Count_sPlot | Count_iNat | |
|---|---|---|---|---|---|
| 0 | -155.0 | 63.0 | Dispersal unit length | 109 | 0 |
| 1 | -155.0 | 69.0 | Dispersal unit length | 66 | 0 |
| 2 | -157.0 | 69.0 | Dispersal unit length | 7 | 2 |
| 3 | -163.0 | 67.0 | Dispersal unit length | 1 | 40 |
| 4 | -157.0 | 65.0 | Dispersal unit length | 1 | 13 |
Load grid cell means at 2 degree resolution:
filename="grid_means_2_deg.csv"
means = pd.read_csv(filename)
means.head()
| x_bin | y_bin | Trait | TraitValue_sPlot | TraitValue_iNat | |
|---|---|---|---|---|---|
| 0 | -157.0 | 69.0 | Dispersal unit length | 0.872919 | 1.898598 |
| 1 | -163.0 | 67.0 | Dispersal unit length | 1.282093 | 1.265103 |
| 2 | -157.0 | 65.0 | Dispersal unit length | 1.267266 | 1.079049 |
| 3 | -153.0 | 65.0 | Dispersal unit length | 1.326229 | 1.485913 |
| 4 | -147.0 | 69.0 | Dispersal unit length | 0.747151 | 0.876784 |
Normalize grid cell averages:
def quantile_norm(df, s1, s2, variables):
# empty data frame to save output:
df_norm = pd.DataFrame()
for v in variables:
# make subset df
sub_exp = df[df['Trait']==v]
sub_exp[s1] = np.exp(sub_exp[s1].copy())
sub_exp[s2] = np.exp(sub_exp[s2].copy())
# determine min and max values
min_quantile = sub_exp[s1].quantile(0.05)
max_quantile = sub_exp[s1].quantile(0.95)
if min_quantile > sub_exp[s2].quantile(0.05):
min_quantile = sub_exp[s2].quantile(0.05)
if max_quantile < sub_exp[s2].quantile(0.95):
max_quantile = sub_exp[s2].quantile(0.95)
sub_exp[s1] = sub_exp[s1].apply(lambda x: (x - min_quantile)/(max_quantile - min_quantile))
sub_exp[s2] = sub_exp[s2].apply(lambda x: (x - min_quantile)/(max_quantile - min_quantile))
df_norm = pd.concat([df_norm, sub_exp])
return df_norm
# normalize original values (exp of ln-values):
pd.options.mode.chained_assignment = None
means = quantile_norm(means, "TraitValue_sPlot", "TraitValue_iNat", trait)
means.head()
| x_bin | y_bin | Trait | TraitValue_sPlot | TraitValue_iNat | Difference | |
|---|---|---|---|---|---|---|
| 0 | -157.0 | 69.0 | Dispersal unit length | 0.133789 | 0.777220 | 0.643432 |
| 1 | -163.0 | 67.0 | Dispersal unit length | 0.315624 | 0.306502 | 0.009122 |
| 2 | -157.0 | 65.0 | Dispersal unit length | 0.307655 | 0.216120 | 0.091535 |
| 3 | -153.0 | 65.0 | Dispersal unit length | 0.340059 | 0.438045 | 0.097986 |
| 4 | -147.0 | 69.0 | Dispersal unit length | 0.091284 | 0.135182 | 0.043898 |
# calculate absolute difference
means['Difference'] = abs(means['TraitValue_iNat'] - means['TraitValue_sPlot'])
means = pd.merge(means, counts, on = ['x_bin', 'y_bin', 'Trait' ])
Remove outliers:
import scipy.stats as stats
z_scores = stats.zscore(means['Difference'])
abs_z_scores = np.abs(z_scores)
filtered_entries = (abs_z_scores < 3)
means = means[filtered_entries]
Plot distribution of the difference between iNautralist and sPlotOpen in each respective cell:
fig,ax = plt.subplots(figsize=(5,5))
plt.hist(means["Difference"], range=[0, 2], color="midnightblue")
plt.title("Absolute difference distribtution")
ax.set(xlabel = "Absolute scaled difference (iNat - sPlotOpen)", ylabel="frequency")
plt.savefig('../Figures/Absolute_difference_distribtution.pdf', bbox_inches='tight')
fig,ax = plt.subplots(figsize=(2,5))
plt.boxplot(means["Count_iNat"])
plt.title("Cell Counts iNaturalist")
ax.set_yscale('log')
plt.savefig('../Figures/Dist_cell_counts.pdf', bbox_inches='tight')
We bin all grid cells by the number of observations they hold. Each bin is defined by the 25% quantiles in the dataset. We can see that the cells in the first bin have only a maximum of 20 observations and exhibit a greater discrepancy from sPlotOpen.
# define 25% quantile bins
means['Counts per grid cell'] = pd.cut(means['Count_iNat'], bins=[0
,means['Count_iNat'].quantile(0.25)
,means['Count_iNat'].quantile(0.50)
,means['Count_iNat'].quantile(0.75)
,means['Count_iNat'].quantile(1)])
fig,ax = plt.subplots(figsize=(7,5))
sns.barplot(x='Counts per grid cell', y='Difference', data=means, ax=ax, color="midnightblue", errcolor = "red")
ax.set(xlabel = "Counts per grid cell in 25% quantiles", ylabel="Absolute difference (iNat - sPlotOpen)")
ax.set(title= "iNat Counts vs. absolute Difference")
plt.savefig('../Figures/diff_vs_counts_all.pdf', bbox_inches='tight')
For each trait individually:
for t in trait:
sub = means[means['Trait']==t]
sub['Counts per grid cell'] = pd.cut(sub['Count_iNat'], bins=[0
,sub['Count_iNat'].quantile(0.25)
,sub['Count_iNat'].quantile(0.50)
,sub['Count_iNat'].quantile(0.75)
,sub['Count_iNat'].quantile(1)])
fig,ax = plt.subplots(figsize=(7,5))
sns.barplot(x='Counts per grid cell', y='Difference', data=sub,
ax=ax, color="midnightblue", errcolor = "red").set(title=t)
ax.set(xlabel = "Counts per grid cell in 25% quantiles", ylabel="Difference (iNat - sPlotOpen)");
Same plot for sPlotOpen plot count
means['Counts per grid cell'] = pd.cut(means['Count_sPlot'], bins=[0
,means['Count_sPlot'].quantile(0.25)
,means['Count_sPlot'].quantile(0.50)
,means['Count_sPlot'].quantile(0.75)
,means['Count_sPlot'].quantile(1)])
fig,ax = plt.subplots(figsize=(7,5))
sns.barplot(x='Counts per grid cell', y='Difference', data=means, ax=ax, color="midnightblue", errcolor = "red")
ax.set(xlabel = "Counts per grid cell in 25% quantiles", ylabel="Difference (iNat - sPlotOpen)")
plt.savefig('../Figures/diff_vs_counts_sPlot.pdf', bbox_inches='tight')
