.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples\apw.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_apw.py: APW Example ====================================== This example demonstrates fitting a penalized spherical spline to the Apparent Polar Wander (APW) path and plotting the result on a world map. .. image:: /_static/thumbnails/apw_thumb.png :class: sphx-glr-thumbimg .. GENERATED FROM PYTHON SOURCE LINES 11-150 .. image-sg:: /auto_examples/images/sphx_glr_apw_001.png :alt: APW Data + Fitted Spherical Spline :srcset: /auto_examples/images/sphx_glr_apw_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Best λ index: 28 Control points (cartesian): [[-0.48799339 -0.18359966 0.85331918] [-0.41541484 -0.15888491 0.89564842] [-0.05778102 0.06411578 0.9962683 ] [-0.07458957 0.22358783 0.97182554] [-0.07418727 0.22529426 0.97146217] [-0.0415261 0.16843219 0.98483815] [-0.02021256 0.17445605 0.98445748] [ 0.06310967 0.32686588 0.94296122] [ 0.02050204 0.37587383 0.92644403] [ 0.00531943 0.39425323 0.91898645] [-0.05242935 0.50251786 0.86297564] [-0.05414767 0.50781569 0.85976232] [-0.08338042 0.52447341 0.84733426] [-0.50457063 0.33962265 0.7937663 ]] | .. code-block:: Python # ---------------------------------------------------- # Imports # ---------------------------------------------------- import numpy as np import pandas as pd import matplotlib.pyplot as plt import geopandas as gpd import spheresmooth as ss from spheresmooth.loads import load_world_map # ---------------------------------------------------- # 1. Load APW spherical data (θ, φ in radians) # ---------------------------------------------------- apw = ss.load_apw() # pandas DataFrame: [t, theta, phi] t = apw.iloc[:, 0].values spherical = apw.iloc[:, 1:3].values # (theta, phi) # ---------------------------------------------------- # 2. Spherical → Cartesian # ---------------------------------------------------- apw_cartesian = ss.spherical_to_cartesian(spherical) # ---------------------------------------------------- # 3. Quantile knots # ---------------------------------------------------- dimension = 15 initial_knots = ss.knots_quantile(t, dimension) lambda_seq = np.exp(np.linspace(np.log(1e-7), np.log(1), 40)) # ---------------------------------------------------- # 4. Penalized spherical spline fit # ---------------------------------------------------- fit = ss.penalized_linear_spherical_spline( t=t, y=apw_cartesian, dimension=dimension, initial_knots=initial_knots, lambdas=lambda_seq ) bic_list = fit["bic_list"] fits = fit["fits"] best_index = np.argmin(bic_list) best_fit = fits[best_index] control_points = best_fit["control_points"] knots = best_fit["knots"] print("Best λ index:", best_index) print("Control points (cartesian):") print(control_points) # ---------------------------------------------------- # 5. Convert control points → (theta, phi) # ---------------------------------------------------- cp_spherical = ss.cartesian_to_spherical(control_points) cp_deg = np.degrees(cp_spherical) cp_df = pd.DataFrame({ "latitude": 90 - cp_deg[:, 0], "longitude": cp_deg[:, 1] }) # ---------------------------------------------------- # 6. Evaluate piecewise geodesic curve # ---------------------------------------------------- ts = np.linspace(0, 1, 2000) curve_cart = ss.piecewise_geodesic(ts, control_points, knots) curve_sph = ss.cartesian_to_spherical(curve_cart) curve_deg = np.degrees(curve_sph) curve_df = pd.DataFrame({ "latitude": 90 - curve_deg[:, 0], "longitude": curve_deg[:, 1] }) # ---------------------------------------------------- # 7. Original APW data → (lat, lon) # ---------------------------------------------------- apw_deg = np.degrees(spherical) apw_df = pd.DataFrame({ "latitude": 90 - apw_deg[:, 0], "longitude": apw_deg[:, 1] }) # ---------------------------------------------------- # 8. Plot world map + APW + spline curve (single figure) # ---------------------------------------------------- fig, ax = plt.subplots(figsize=(12, 10)) # World map world_path = load_world_map() world = gpd.read_file(world_path) world.plot(ax=ax, color="antiquewhite", edgecolor="grey") # Raw APW data ax.scatter( apw_df["longitude"], apw_df["latitude"], s=5, color="black", label="APW data" ) # Control points ax.scatter( cp_df["longitude"], cp_df["latitude"], s=60, color="blue", marker="s", label="Control points" ) # Fitted spline curve ax.plot( curve_df["longitude"], curve_df["latitude"], c="red", linewidth=2, label="Spline curve" ) # --------------------------- # Zoom region # --------------------------- ax.set_xlim(60, 230) ax.set_ylim(45, 95) ax.set_xlabel("Longitude") ax.set_ylabel("Latitude") ax.set_title("APW Data + Fitted Spherical Spline") ax.legend() plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 47.499 seconds) .. _sphx_glr_download_auto_examples_apw.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: apw.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: apw.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: apw.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_