finishing f4, f6, fs2

README.md

@@ -54,7 +54,7 @@ The figures are in [`article_figs`](article_figs/). The figures with their respe
 
 **Supplementary Information Figures**
 - **Figure S1**: [`supplementary_material_plots.ipynb`](analysis/supplementary_material_plots.ipynb)
-- **Figure S2**: [`vertical_histograms.py`](analysis/vertical_histograms.py)
+- **Figure S2**: [`Map_origin_particles_surface.ipynb`](Map_origin_particles_surface.ipynb)
 
 ## Data Availability
 The data will be made available with its own DOI, when the review process is over.

analysis/analysis_functions.py

@@ -156,7 +156,7 @@ def ridge_plot(data, xlabel ='', title='', bins=128, h_space=-0.5, alpha=1,
 cmap_bathy = sns.mpl_palette('cmo.ice', n_colors=8, as_cmap=True)
 
 
-def bathymetry_plot(figsize=(13, 7),alpha=1., ):
+def bathymetry_plot(figsize=(13, 7),alpha=1., land_zorder=5):
 
     fig = plt.figure(figsize=figsize)
     ax = plt.axes(projection=ccrs.PlateCarree())
@@ -172,7 +172,8 @@ def bathymetry_plot(figsize=(13, 7),alpha=1., ):
     ax.add_feature(cartopy.feature.NaturalEarthFeature(category='physical',
                                                        name='land', 
                                                        scale='110m'),
-                                                       color='black')
+                                                       color='black',
+                                                       zorder=land_zorder)
 
     gl = ax.gridlines(draw_labels=True)
     gl.right_labels = False

analysis/size_distribution.py

@@ -36,6 +36,7 @@
     # print('Computing fragmentation timescale: ', ft)
     local_path = f'/storage/shared/oceanparcels/output_data/data_Claudio/hc13_3/hc13_{ft}.zarr'
     sim = xr.open_zarr(local_path)
+    sim = sim.where(sim.time >= np.datetime64('2007-01-01'), drop=True) # analysis stops at 2007-01-01
 
     # loading fields as np arrays to do some nasty indexing
     radiuses = sim['radius'].values
@@ -167,75 +168,13 @@ def inverse(x):
 fig.savefig('../article_figs/ECDF_surface', dpi=300,
             facecolor=(1, 0, 0, 0))
 
-
-# %% Results ECDF
-frag_into_NPs = np.load('../data/frag_into_NPs.npy', allow_pickle=True)[()]
-
-fig, ax = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
-
-ax[0].axvline(initial_depth, ls=':', color='k')
-ax[0].text(initial_depth - 200, 0.65, r'Sampling Depth', fontsize=6, color='k', rotation=-90)
-ax[0].axvline(0, ls=':', color='k')
-ax[0].text(0-200, 0.65, r'Surface', fontsize=6, color='k', rotation=-90)
-
-ax[1].axvline(1e-6, ls=':', color='black')
-ax[1].axvline(1e-4, ls=':', label=r"Fragmentation limit", color='red')
-ax[1].text(1e-6, 0.08, r"1 $\mu m$ Limit", fontsize=6, color='k', rotation=-90)
-ax[1].text(1.1e-4, 0.01, r"Fragmentation Limit", fontsize=6, color='r', rotation=-90)
-
-for j, ft in enumerate(simulations[::-1]):
-
-    x, y = funk.ecdf(frag_into_NPs[ft]['depths'], normalized=True,
-                     invert=False)
-    ax[0].plot(x, y, drawstyle='steps-post', label=f'$\lambda_f$ = {ft} days')
-
-    x, y = funk.ecdf(frag_into_NPs[ft]['particle_index'], normalized=True,
-                     invert=False)
-
-    x, y = funk.ecdf(surface_events[ft]['radius'], normalized=True)
-    ax[1].plot(x, y, drawstyle='steps-post')
-
-
-handles, labels = ax[0].get_legend_handles_labels()
-handles = handles[::-1]
-labels = labels[::-1]
-
-ax[0].legend(handles, labels, fontsize=7, shadow=True, ncol=2,
-         loc='best')
-
-ax[1].legend(handles, labels, fontsize=7, shadow=True, ncol=1,
-         loc='best')
-
-ax[0].set_xlabel('$R < 1\ \mu m$ Fragmentation Depth, $z$ [m]')
-ax[0].set_ylabel(r'ECDF: $P(x \leq z)$')
-
-ax[1].semilogx()
-ax[1].set_xlabel('Surface Particles Radius, $R$ [m]')
-ax[1].set_ylabel(r'ECDF: $P(x \leq R)$')
-
-gridy = np.linspace(0, 1, 11)
-gridx = [500, 1000] + [i for i in range(2000, 10000, 2000)]
-
-ax[0].set_yticks(gridy)
-ax[1].set_yticks(gridy)
-
-ax[0].grid()
-ax[1].grid()
-
-# ax[0].text(-5500, 0.98, r'A', fontsize=12,
-#                ha='right')
-# ax[1].text(1e-7, 0.98, r'B', fontsize=12,
-#                ha='right')
-
-# fig.savefig('../article_figs/ECDFs2.png', dpi=300,
-#             facecolor=(1, 0, 0, 0))
 # %% FIGURE 6 -
 fig, ax = plt.subplots(1, 1, figsize=(4, 3.5), tight_layout=True)
 
-ax.axvline(initial_depth, ls=':', color='k')
+ax.axvline(initial_depth, ls='-', color='k', lw=1)
 ax.text(initial_depth - 200, 0.65, r'Sampling Depth', fontsize=6, color='k', rotation=-90)
-ax.axvline(0, ls=':', color='k')
-ax.text(0-200, 0.65, r'Surface', fontsize=6, color='k', rotation=-90)
+ax.axvline(0, ls='-', color='k', lw=1)
+ax.text(0-220, 0.7, r'Surface', fontsize=6, color='k', rotation=-90)
 
 colors = plt.get_cmap('tab10').colors
 line_styles = ['-', '--', '-.', ':']
@@ -259,7 +198,7 @@ def inverse(x):
 
 gridy = np.linspace(0, 1, 11)
 ax.set_yticks(gridy)
-ax.grid()
+ax.grid(linestyle=':')
 
 fig.savefig('../article_figs/Figure6.png', dpi=300,
             facecolor=(1, 0, 0, 0))
@@ -268,10 +207,10 @@ def inverse(x):
 
 fig, ax = plt.subplots(1, 1, figsize=(4, 3.5), tight_layout=True)
 
-ax.axvline(1e-6, ls=':', color='black')
-ax.axvline(1e-4, ls=':', color='red')
+ax.axvline(1e-6, ls='-', lw=1, color='black')
+ax.axvline(1e-4, ls='-', lw=1, color='black')
 ax.text(1e-6, 0.08, r"1 $\mu m$ Limit", fontsize=6, color='k', rotation=-90)
-ax.text(1.1e-4, 0.01, r"Fragmentation Limit", fontsize=6, color='r', rotation=-90)
+ax.text(1.1e-4, 0.01, r"Fragmentation Limit", fontsize=6, color='black', rotation=-90)
 
 colors = plt.get_cmap('tab10').colors
 line_styles = ['-', '--', '-.', ':']
@@ -295,8 +234,10 @@ def inverse(x):
 
 gridy = np.linspace(0, 1, 11)
 ax.set_yticks(gridy)
-ax.grid()
+ax.grid(linestyle=':')
 
 fig.savefig('../article_figs/Figure4.png', dpi=300,
             facecolor=(1, 0, 0, 0))
 
+
+# %%

analysis/vertical_histograms.py

@@ -15,7 +15,7 @@
 from datetime import datetime
 import analysis_functions as funk
 
-run_for_loop = True
+run_for_loop = False
 
 # Define initial conditions
 initial_depth = -5000  # int(sys.argv[1])  # 5 # 60 # 5179
@@ -158,8 +158,7 @@
                ha='right')
     ax[j].set_yticks([-5500, -2500, 0])
     ax[j].grid()
-
-ax[1].set_ylabel('Depth (m)')
+    ax[j].set_ylabel('Depth (m)')
 
 fig.colorbar(im, ax=ax[-1], orientation='horizontal',
              extend='max', label='Depth Probability of Nanoplastics')
@@ -172,7 +171,8 @@
 
 # %% Depth vs displacement Plot
 # '-', '--', '-.', ':', 'None', ' ', '', 'solid', 'dashed', 'dashdot', 'dotted'
-marker = itertools.cycle(('v', 'h', 'd', 'o', 'X', 'P', '^', 's'))
+marker = itertools.cycle(('o', 'd', 'X', 's'))
+# marker = itertools.cycle(('v', 'h', 'd', 'o', 'X', 'P', '^', 's'))
 
 fig = plt.figure(figsize=(8, 4))
 gs = GridSpec(2, 2, width_ratios=[4, 1], height_ratios=[1, 4], wspace=0.01,
@@ -306,48 +306,4 @@ def inverse(x):
 fig.savefig('../article_figs/ECDF_nanoparticles', dpi=300,
             facecolor=(1, 0, 0, 0))
 
-# %% FIGURE S2 - Maps of fragmenting location
-marker = itertools.cycle(('v', 'h', 'd', 'o', 'X', 'P', '^', 's'))
-fig, ax = funk.bathymetry_plot(alpha=0.1)
-
-for j, ft in enumerate(simulations[::-1]):
-    ax.scatter(frag_into_NPs[ft]['lon'], frag_into_NPs[ft]['lat'], zorder=2,
-               s=20,
-               label=f"$\lambda_f$ = {ft} days",
-               marker=next(marker))
-
-ax.scatter(origin[0], origin[1], zorder=5,
-           label='Sampling Location', marker='*', s=80, edgecolors='k')
-
-for r in range(1,10):
-    circle_points = geodesic.Geodesic().circle(lon=origin[0], lat=origin[1],
-                                               radius=r*1e6,
-                                               n_samples=100,
-                                               endpoint=False)
-    geom = shapely.geometry.Polygon(circle_points)
-    ax.add_geometries((geom,), crs=ccrs.PlateCarree(), facecolor='none',
-                      edgecolor='black', linewidth=1., zorder=3, ls='--')
-
-# ax.set_title('Where do the particles fragment into Nanoparticles?')
-handles, labels = ax.get_legend_handles_labels()
-handles = handles[::-1]
-labels = labels[::-1]
-
-ax.text(3.3, -24.5, r"$1,000$ km", fontsize=5)
-ax.text(-14.5, -30.5, r"$2,000$ km", fontsize=5, rotation=70)
-ax.text(-25., -30.5, r"$3,000$ km", fontsize=5, rotation=70)
-ax.text(-35.6, -30.5, r"$4,000$ km", fontsize=5, rotation=70)
-ax.text(-46.6, -31.5, r"$5,000$ km", fontsize=5, rotation=70)
-ax.text(-63.5, -48.5, r"$6,000$ km", fontsize=5, rotation=75)
-ax.text(-79.3, -52.5, r"$7,000$ km", fontsize=5, rotation=62)
-ax.text(-89.5, -45.3, r"$8,000$ km", fontsize=5, rotation=53)
-ax.text(-89.5, -27.3, r"$9,000$ km", fontsize=5, rotation=55)
-
-ax.set_extent([-97, 60, -65, 0], crs=ccrs.PlateCarree())
-
-ax.legend(handles, labels, ncols=3, fontsize=9, shadow=True, loc='upper left')
-
-fig.savefig('../article_figs/FigureS2.png', dpi=300,
-            facecolor=(1, 0, 0, 0))
-
 # %%