AntSimi / py-eddy-tracker

Eddy identification and tracking
https://py-eddy-tracker.readthedocs.io/en/latest/
GNU General Public License v3.0
123 stars 52 forks source link

Help in generating representative contour and save as .nc #248

Open hafez-ahmad opened 2 months ago

hafez-ahmad commented 2 months ago

Thank you so much for creating wonderfull Package. I am trying to reproduce following

https://py-eddy-tracker.readthedocs.io/en/latest/python_module/10_tracking_diagnostics/pet_pixel_used.html#sphx-glr-python-module-10-tracking-diagnostics-pet-pixel-used-py

Could you please help me on how to add some representative contour lines on top this plot? and how I can save as nc/ geotif file so I could further work on this.

Here my code

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import LogNorm

step = 0.125
bins = ((75, 100, step), (2, 23, step))
kwargs_pcolormesh = dict(
    cmap="terrain_r", vmin=0, vmax=0.75, factor=1 / a.nb_days, name="count"
)

fig, ((ax_a, ax_c), (ax_all, ax_ratio)) = plt.subplots(2, 2, figsize=(12, 12), subplot_kw={'projection': ccrs.PlateCarree()})

ax_a.set_title("Anticyclonic frequency")
ax_c.set_title("Cyclonic frequency")
ax_all.set_title("All eddies frequency")
ax_ratio.set_title("Ratio Cyclonic / Anticyclonic")

# Add labels
ax_a.text(-0.1, 1.1, 'A', transform=ax_a.transAxes, fontsize=16, fontweight='bold', va='top')
ax_c.text(-0.1, 1.1, 'B', transform=ax_c.transAxes, fontsize=16, fontweight='bold', va='top')
ax_all.text(-0.1, 1.1, 'C', transform=ax_all.transAxes, fontsize=16, fontweight='bold', va='top')
ax_ratio.text(-0.1, 1.1, 'D', transform=ax_ratio.transAxes, fontsize=16, fontweight='bold', va='top')

# Function to setup map features
def setup_map(ax):
    ax.set_extent([75, 100, 2, 23], crs=ccrs.PlateCarree())
    ax.coastlines()
    ax.add_feature(cfeature.LAND, color='gray')
    gl = ax.gridlines(draw_labels=True)
    gl.top_labels = False
    gl.right_labels = ax in [ax_c, ax_ratio]  # Show right labels only on the right plots
    gl.left_labels = ax in [ax_a, ax_all]  # Show left labels only on the left plots
    gl.bottom_labels = True

# Setup map features for each subplot
for ax in (ax_a, ax_c, ax_all, ax_ratio):
    setup_map(ax)

# Count pixel used for each contour
g_a = a.grid_count(bins, intern=True)
m_a = g_a.display(ax_a, **kwargs_pcolormesh)
g_c = c.grid_count(bins, intern=True)
m_c = g_c.display(ax_c, **kwargs_pcolormesh)

# Add a common color bar for the first two plots (ax_a and ax_c)
#cbar_ax = fig.add_axes([0.2, 0.53, 0.6, 0.02])  # Adjust the position and size of the color bar
#fig.colorbar(m_a, cax=cbar_ax, orientation='horizontal')

# Compute a ratio Cyclonic / Anticyclonic
ratio = g_c.vars["count"] / g_a.vars["count"]

# Mask manipulation to be able to sum the 2 grids
m_c = g_c.vars["count"].mask
m = m_c & g_a.vars["count"].mask
g_c.vars["count"][m_c] = 0
g_c.vars["count"] += g_a.vars["count"]
g_c.vars["count"].mask = m

m_all = g_c.display(ax_all, **kwargs_pcolormesh)
cbar_ax_all = fig.add_axes([0.92, 0.55, 0.02, 0.35])  # Adjust the position and size of the color bar for ax_all
fig.colorbar(m_all, cax=cbar_ax_all, orientation='vertical')

g_c.vars["count"] = ratio
m_ratio = g_c.display(
    ax_ratio, name="count", norm=LogNorm(vmin=0.1, vmax=10), cmap="coolwarm_r"
)
cbar_ax_ratio = fig.add_axes([0.92, 0.13, 0.02, 0.34])  # Adjust the position and size of the color bar for ax_ratio
fig.colorbar(m_ratio, cax=cbar_ax_ratio, orientation='vertical')

plt.tight_layout(rect=[0, 0, 0.9, 1])  # Adjust the layout to make space for the color bars

# Save figure
plt.savefig('Frequency.jpg', dpi=500)
AntSimi commented 2 months ago

To my knowledge, matplotlib does not allow you to create geotiffs.

Could you please help me on how to add some representative contour lines on top this plot? Contour of which data?

hafez-ahmad commented 2 months ago

Thank you so much for your response. @AntSimi for all four figures of Anticyclonic frequency, Cyclonic frequency All eddies frequency, Ratio Cyclonic / Anticyclonic. I like to add contour on top of frequency.

AntSimi commented 2 months ago

Take a look at this method RegularGridDataset.contour

hafez-ahmad commented 2 months ago

Could you please write demo like you created deom on frequency ?

Thank you

AntSimi commented 2 months ago

Use is similar to grid display method, keyword argument are listed on documentation and on matplotlib documentation linked in pet contour documentation. Use could look like this if you want contour for some levels:

replace_by_my_grid.contour(replace_by_my_ax, 'replace_by_my_var', levels=[-3500,-2500,-1500])

Use could look like this if you want contour for a number of levels:

replace_by_my_grid.contour(replace_by_my_ax, 'replace_by_my_var', levels=5)
hafez-ahmad commented 1 week ago

Thank you so much. Could you please show any way to calculate kinetic energy and baroclinic and barotropic instability ?