AntSimi / py-eddy-tracker

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

The question about the data #214

Open SomnusQue opened 10 months ago

SomnusQue commented 10 months ago
SomnusQue commented 10 months ago

Dear authors, I am wondering the data structure about this project. Now, I got a SLA.nc file which could not run the pet_sla_and_adt.py, the reason why that was Coordinates in RegularGridDataset must be strictly increasing. I have been sort the longitude and latitude about this data, should I make sure the data is arithmetic progression? Best wishes.

AntSimi commented 10 months ago

could you share ncdump -h of your file?

SomnusQue commented 10 months ago

Of course, this is our file. Now, the new error that was Cannot convert from 'dimensionless' (dimensionless) to 'meter' ([length])..I'm sorry to trouble you again. Would you please see what the problem is? By the way, our file don't contains adt data. We just want to detect eddy by SLA.. Best wishes pred_mra5_20190101_zos.nc.zip

SomnusQue commented 10 months ago

This is the code which we changed.

"""
Eddy detection on SLA
=============================

"""
from datetime import datetime

from matplotlib import pyplot as plt

from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset

def start_axes(title):
    fig = plt.figure(figsize=(13, 5))
    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
    # ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
    ax.set_aspect("equal")
    ax.set_title(title)
    return ax

def update_axes(ax, mappable=None):
    ax.grid()
    if mappable:
        plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))

g = RegularGridDataset(
    data.get_demo_path("/Users/Maybe/Desktop/pred_mra5_20190101_zos.nc"),
    "longitude",
    "latitude",
)

g.add_uv("prediction", "ugosa", "vgosa")
wavelength = 400
print('here')
g.copy("prediction", "sla_raw")

g.bessel_high_filter("prediction", wavelength)
date = datetime(2019, 1, 1)

kwargs_a_sla = dict(
    lw=0.5, label="Anticyclonic SLA ({nb_obs} eddies)", ref=-10, color="g"
)
kwargs_c_sla = dict(lw=0.5, label="Cyclonic SLA ({nb_obs} eddies)", ref=-10, color="b")

a_sla, c_sla = g.eddy_identification("prediction", "ugosa", "vgosa", date, 0.002)

ax = start_axes(f"SLA (m) filtered ({wavelength}km)" + str(date.date()))
m = g.display(ax, "sla", vmin=-0.15, vmax=0.15)
a_sla.display(ax, **kwargs_a_sla), c_sla.display(ax, **kwargs_c_sla)
ax.legend(), update_axes(ax, m)

ax = start_axes("SLA (m)")
m = g.display(ax, "sla_raw", vmin=-0.15, vmax=0.15)
a_sla.display(ax, **kwargs_a_sla), c_sla.display(ax, **kwargs_c_sla)
ax.legend(), update_axes(ax, m)

ax = start_axes("Eddies detected")
a_sla.display(ax, **kwargs_a_sla)
c_sla.display(ax, **kwargs_c_sla)
ax.legend()
update_axes(ax)
plt.show()
AntSimi commented 9 months ago

First data.get_demo_path method is only to get demonstration sample, you don't need it for your own data. You need to specify height and speed unit in eddy_identification method if there is no attribute unit in your netcdf.