Example of sampling

Example of sampling#

Here we illustrate an example of how to generate samples using our small Python package smpsite. Many of the sampling capabilities are directly imported from pmagpy.

import pmagpy.pmag as pmag
import pmagpy.ipmag as ipmag
import matplotlib.pyplot as plt
import smpsite as smp

import numpy as np
np.random.seed(666)

%matplotlib inline

Define the parameters of the sample#

We define a parameter class defined inside smpsite that allow us to specify all the parameters of interest at the same time. These parameters include:

N: Total number of sites
n0: Number of samples per site
kappa_within_site: Concentration parameter of the Fisher distribution within each site
site_lat: Latitude of the site
site_long: Longitude of the site
outlier_rate: Fraction of samples that are randomly taken from the uniform distribution in the sphere, with 0.0 representing no outliers and 1.0 representing all outliers.
secular_method: Method used to sample the VGPs. It includes G, tk03 and Fisher.
kappa_secular: In case of secular_method="Fisher", it specifies the concentration parameter of the Fisher distribution.

params0 = smp.Params(N=10,
                     n0=5,
                     kappa_within_site=100,
                     site_lat=10, 
                     site_long=0,
                     outlier_rate=0.10,
                     secular_method="G",
                     kappa_secular=None)

Generate the sample#

Once the parameters are specified, we can simulate the sample data:

%%time

df_sample = smp.generate_samples(params0)
df_sample.head()

CPU times: user 16.6 ms, sys: 1.05 ms, total: 17.6 ms
Wall time: 17.1 ms

	vgp_long	vgp_lat	vgp_dec	vgp_inc	is_outlier
0	126.657018	82.567616	5.984268	10.932745	0
1	127.014390	76.856329	10.467458	3.924856	0
2	91.215608	78.611284	11.548712	18.618003	0
3	103.137241	78.865732	10.929083	14.380461	0
4	68.461601	-57.543217	150.016322	5.433526	1

Plot the simulated sample data#

We can plot the declination and inclination of all of the simulated samples.

plt.figure()
ipmag.plot_net(1)
ipmag.plot_di(dec=df_sample.vgp_dec.values, inc=df_sample.vgp_inc.values)
plt.title('declination/inclination of samples')
plt.show()

../_images/47f57c4fdec3c73f8316ad26af0f4288f5aaef92980f0f61630479a004fcb2ab.png

The declination/inclination can be plotted as colored them by site

# Create a color map based on unique 'sample_site' values
unique_sites = df_sample['sample_site'].unique()
colors = plt.cm.viridis(np.linspace(0, 1, len(unique_sites)))
color_map = dict(zip(unique_sites, colors))

ipmag.plot_net(1)

# Use a set to track sites that have already been added to the legend
added_legend_sites = set()

# Iterate through the dataframe and plot each value with its unique color
for index, row in df_sample.iterrows():
    site = row['sample_site']
    color = color_map[site]
    # Using your given column names for declination and inclination
    ipmag.plot_di(row['vgp_dec'], row['vgp_inc'], color=color)
    
    # If the site hasn't been added to the legend, add it
    if site not in added_legend_sites:
        # Use marker='o' for circles, linestyle='' to remove line, and convert site to int for the label
        plt.plot([], [], 'o', label=int(site), color=color, linestyle='')
        added_legend_sites.add(site)

# Show the legend
plt.legend(title='sites',loc='best')
plt.title('declination/inclination of samples\ncolored by site')
plt.show()

../_images/238ae45e6d7758e7754774d75d0381ce521aea5bdeaeb6d1afebb4a67021d0e9.png

Calculate pole position and calculate angular error#

mean_pole = smp.estimate_pole(df_sample, params0, ignore_outliers="True")
print(mean_pole)
angular_error = pmag.angle([mean_pole['pole_dec'], mean_pole['pole_inc']], [0, 90])
print('The angular error is: ', round(angular_error[0],1), 'degrees')

{'pole_dec': np.float64(112.27111493549478), 'pole_inc': np.float64(87.17580406361836), 'S2_vgp': np.float64(256.9146243330398), 'total_samples': np.float64(48.0), 'samples_per_site': 5, 'alpha95': np.float64(9.814841392668617)}
The angular error is:  2.8 degrees

/home/runner/work/PaleoSampling/PaleoSampling/smpsite/smpsite/estimate.py:119: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
  df_site = df.groupby('sample_site').apply(lambda row : pd.Series(robust_fisher_mean(row.vgp_dec.values, row.vgp_inc.values)))

Plot the pole position#

map_axis = ipmag.make_orthographic_map(central_longitude=0,central_latitude=90)
ipmag.plot_pole(map_axis, mean_pole['pole_dec'], mean_pole['pole_inc'], mean_pole['alpha95'], color='r', marker='s', label='Mean Pole')

/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:6195: DeprecationWarning: invalid escape sequence \d
  pattern = re.compile('[4][-]\d')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:6202: DeprecationWarning: invalid escape sequence \d
  pattern = re.compile('[7][-]\d')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7886: DeprecationWarning: invalid escape sequence \h
  sf.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7891: DeprecationWarning: invalid escape sequence \h
  f.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7896: DeprecationWarning: invalid escape sequence \h
  fI.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7902: DeprecationWarning: invalid escape sequence \h
  cr.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7908: DeprecationWarning: invalid escape sequence \h
  fsp.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7942: DeprecationWarning: invalid escape sequence \h
  f.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7943: DeprecationWarning: invalid escape sequence \h
  sf.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7944: DeprecationWarning: invalid escape sequence \h
  fI.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7946: DeprecationWarning: invalid escape sequence \h
  cr.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:7948: DeprecationWarning: invalid escape sequence \h
  fsp.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8060: DeprecationWarning: invalid escape sequence \h
  f.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8061: DeprecationWarning: invalid escape sequence \h
  sf.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8062: DeprecationWarning: invalid escape sequence \h
  fI.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8063: DeprecationWarning: invalid escape sequence \e
  f.write('\end{longtable}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8064: DeprecationWarning: invalid escape sequence \e
  sf.write('\end{longtable}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8065: DeprecationWarning: invalid escape sequence \e
  fI.write('\end{longtable}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8066: DeprecationWarning: invalid escape sequence \e
  f.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8067: DeprecationWarning: invalid escape sequence \e
  sf.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8068: DeprecationWarning: invalid escape sequence \e
  fI.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8070: DeprecationWarning: invalid escape sequence \h
  fsp.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8071: DeprecationWarning: invalid escape sequence \e
  fsp.write('\end{longtable}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8072: DeprecationWarning: invalid escape sequence \e
  fsp.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8074: DeprecationWarning: invalid escape sequence \h
  cr.write('\hline\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8075: DeprecationWarning: invalid escape sequence \e
  cr.write('\end{longtable}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:8076: DeprecationWarning: invalid escape sequence \e
  cr.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:11626: DeprecationWarning: invalid escape sequence \w
  'LP-*[\w]*-ARM')]
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:11633: DeprecationWarning: invalid escape sequence \s
  cond = spec_df.method_codes.str.contains('(^|[\s\:])LT-PTRM')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12483: DeprecationWarning: invalid escape sequence \d
  directions_out.write('\documentclass{article}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12489: DeprecationWarning: invalid escape sequence \e
  directions_out.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12504: DeprecationWarning: invalid escape sequence \d
  intensities_out.write('\documentclass{article}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12510: DeprecationWarning: invalid escape sequence \e
  intensities_out.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12542: DeprecationWarning: invalid escape sequence \d
  info_out.write('\documentclass{article}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12548: DeprecationWarning: invalid escape sequence \e
  info_out.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12605: DeprecationWarning: invalid escape sequence \d
  info_out.write('\documentclass{article}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12617: DeprecationWarning: invalid escape sequence \e
  info_out.write('\end{landscape}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12618: DeprecationWarning: invalid escape sequence \e
  info_out.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12681: DeprecationWarning: invalid escape sequence \d
  info_out.write('\documentclass{article}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:12689: DeprecationWarning: invalid escape sequence \e
  info_out.write('\end{document}\n')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:13658: DeprecationWarning: invalid escape sequence \c
  plt.ylabel('$\chi$ ($\mu$SI)')
/opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:13678: DeprecationWarning: invalid escape sequence \c
  plt.ylabel('$\chi$ ($\mu$SI)')

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[7], line 1
----> 1 map_axis = ipmag.make_orthographic_map(central_longitude=0,central_latitude=90)
      2 ipmag.plot_pole(map_axis, mean_pole['pole_dec'], mean_pole['pole_inc'], mean_pole['alpha95'], color='r', marker='s', label='Mean Pole')

File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/pmagpy/ipmag.py:1868, in make_orthographic_map(central_longitude, central_latitude, figsize, add_land, land_color, land_edge_color, add_ocean, ocean_color, grid_lines, lat_grid, lon_grid)
   1866     return
   1867 fig = plt.figure(figsize=figsize)
-> 1868 map_projection = LowerThresholdOrthographic(
   1869     central_longitude=central_longitude, central_latitude=central_latitude)
   1870 ax = plt.axes(projection=map_projection)
   1871 ax.set_global()

File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/cartopy/crs.py:2069, in Orthographic.__init__(self, central_longitude, central_latitude, globe)
   2067 self._x_limits = mins[0], maxs[0]
   2068 self._y_limits = mins[1], maxs[1]
-> 2069 self.threshold = np.diff(self._x_limits)[0] * 0.02

AttributeError: can't set attribute

<Figure size 800x800 with 0 Axes>