search

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

eddy identification issue #199

Closed MaYuMMMMM closed 1 year ago

MaYuMMMMM commented 1 year ago

Code for reproduction

from py_eddy_tracker import start_logger
from matplotlib import pyplot as plt
from py_eddy_tracker.dataset.grid import RegularGridDataset
from datetime import datetime

start_logger().setLevel("DEBUG")  # Available options: ERROR, WARNING, INFO, DEBUG

grid_name, lon_name, lat_name = (
    "xinnrt.nc",
    "longitude",
    "latitude",
)

h = RegularGridDataset(grid_name, lon_name, lat_name)

h.bessel_high_filter("adt", 700, order=1)

date = datetime(2023, 4, 16)

a, c = h.eddy_identification(
    "adt",
    "ugos",
    "vgos",  # Variables used for identification
    date,  # Date of identification
    0.002,  # step between two isolines of detection (m)
    pixel_limit=(5, 2000),  # Min and max pixel count for valid contour
    shape_error=55,  # Error max (%) between ratio of circle fit and contour
)

fig = plt.figure(figsize=(15, 7))
ax = fig.add_axes([0.03, 0.03, 0.94, 0.94])
ax.set_title("Eddies detected -- Cyclonic(red) and Anticyclonic(blue)")
ax.set_ylim(15, 25)
ax.set_xlim(-115, -105)
ax.set_aspect("equal")
a.display(ax, color="b", linewidth=0.5)
c.display(ax, color="r", linewidth=0.5)
ax.grid()
fig.savefig("eddies.png")

Actual outcome

WARNING 2023-04-17 18:12:52,753 grid.__init__ :
        We assume pixel position of grid is centered for xinnrt.nc
        DEBUG 2023-04-17 18:12:52,753   grid.   load_general_features :
                Load general feature from xinnrt.nc
        DEBUG 2023-04-17 18:12:52,763   grid.   bessel_high_filter :      
                Run filtering with wavelength of 700 km and order of 1 ...
        DEBUG 2023-04-17 18:12:52,764   grid.   grid :
                Load adt from xinnrt.nc
Remain  0:00:00.001246 ETA  2023-04-17 18:12:52.834959 current kernel size : (53, 51) Step : 12/13    
        DEBUG 2023-04-17 18:12:52,835   grid.   bessel_high_filter :
                Filtering done
        DEBUG 2023-04-17 18:12:52,836   grid.   grid :
                Load ugos from xinnrt.nc
        DEBUG 2023-04-17 18:12:52,841   grid.   grid :
                Load vgos from xinnrt.nc
INFO 2023-04-17 18:12:53,046 grid.eddy_identification :
        We will apply on step a factor to be coherent with grid : 1.000000
        DEBUG 2023-04-17 18:12:53,048   grid.   eddy_identification :
                Levels from -0.169624 to 0.076177
INFO 2023-04-17 18:12:53,048 eddy_feature.__init__ :
        Start computing iso lines
        DEBUG 2023-04-17 18:12:53,078   eddy_feature.   __init__ :
                X shape : (33,)
        DEBUG 2023-04-17 18:12:53,080   eddy_feature.   __init__ :
                Y shape : (13,)
        DEBUG 2023-04-17 18:12:53,081   eddy_feature.   __init__ :
                Z shape : (33, 13)
INFO 2023-04-17 18:12:53,082 eddy_feature.__init__ :
        Start computing iso lines with 125 levels from -0.170000 to 0.078000 ...
INFO 2023-04-17 18:12:53,180 eddy_feature.__init__ :
        Finish computing iso lines
INFO 2023-04-17 18:12:53,183 eddy_feature.__init__ :
        Repair 2 closed contours and 0 almost closed contours / 52 contours
        DEBUG 2023-04-17 18:12:53,185   grid.   eddy_identification :
                doing collection 68, contour value -0.0340, 1 paths
        DEBUG 2023-04-17 18:12:53,238   grid.   eddy_identification :
                doing collection 69, contour value -0.0320, 1 paths
        DEBUG 2023-04-17 18:12:53,239   grid.   eddy_identification :
                doing collection 70, contour value -0.0300, 1 paths
        DEBUG 2023-04-17 18:12:53,316   grid.   eddy_identification :
                doing collection 71, contour value -0.0280, 1 paths
        DEBUG 2023-04-17 18:12:53,317   grid.   eddy_identification :
                doing collection 72, contour value -0.0260, 2 paths
        DEBUG 2023-04-17 18:12:53,322   grid.   eddy_identification :
                doing collection 73, contour value -0.0240, 1 paths
        DEBUG 2023-04-17 18:12:53,323   grid.   eddy_identification :
                doing collection 74, contour value -0.0220, 1 paths
        DEBUG 2023-04-17 18:12:53,326   grid.   eddy_identification :
                doing collection 95, contour value 0.0200, 1 paths
Traceback (most recent call last):
  File "d:/海冰数据/面向全时空/涡旋演变/nc数据/woxuan/woxuan.py", line 22, in <module>
    a, c = h.eddy_identification(
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 849, in eddy_identification
    ) = self.get_uavg(
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 996, in get_uavg
    max_average_speed = self.speed_coef_mean(original_contour)
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 1887, in speed_coef_mean
    return mean_on_regular_contour(
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\numba\core\dispatcher.py", line 468, in _compile_for_args
    error_rewrite(e, 'typing')
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\numba\core\dispatcher.py", line 409, in error_rewrite
    raise e.with_traceback(None)
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Failed in nopython mode pipeline (step: nopython frontend)
Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<built-in function getitem>) found for signature:

 >>> getitem(bool, UniTuple(int64 x 2))

There are 22 candidate implementations:
      - Of which 22 did not match due to:
      Overload of function 'getitem': File: <numerous>: Line N/A.
        With argument(s): '(bool, UniTuple(int64 x 2))':
       No match.

During: typing of intrinsic-call at C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\generic.py (280)

File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\generic.py", line 280:
def interp2d_bilinear(x_g, y_g, z_g, m_g, x, y):
    <source elided>
            if not masked:
                if m_g[i0, j0] or m_g[i0, j1] or m_g[i1, j0] or m_g[i1, j1]:
                ^

During: resolving callee type: type(CPUDispatcher(<function interp2d_bilinear at 0x000001D62323CA60>))
During: typing of call at C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\generic.py (201)

During: resolving callee type: type(CPUDispatcher(<function interp2d_bilinear at 0x000001D62323CA60>))
During: typing of call at C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\generic.py (201)

File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\generic.py", line 201:
def interp2d_geo(x_g, y_g, z_g, m_g, x, y, nearest=False):
    <source elided>
    else:
        return interp2d_bilinear(x_g, y_g, z_g, m_g, x, y)
        ^

During: resolving callee type: type(CPUDispatcher(<function interp2d_geo at 0x000001D62323C700>))
During: typing of call at C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\dataset\grid.py (129)

During: resolving callee type: type(CPUDispatcher(<function interp2d_geo at 0x000001D62323C700>))
During: typing of call at C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\dataset\grid.py (129)

File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 129:
def mean_on_regular_contour(
    <source elided>
    x_new, y_new = uniform_resample(x_val, y_val, num_fac, `fixed_size)`
    values = interp2d_geo(x_g, y_g, z_g, m_g, x_new[1:], y_new[1:])
    ^

image

data.zip

AntSimi commented 1 year ago

Which version of python/numpy/numba did you use?

MaYuMMMMM commented 1 year ago

Python 3.8.3 / numpy 1.20.0 /numba 0.55.2

MaYuMMMMM commented 1 year ago

Code for reproduction

from matplotlib import pyplot as plt
from numpy import arange, outer
import py_eddy_tracker_sample
from py_eddy_tracker.data import get_demo_path
from py_eddy_tracker.observations.network import NetworkObservations
from py_eddy_tracker.observations.observation import EddiesObservations, Table
from py_eddy_tracker.observations.tracking import TrackEddiesObservations

kwargs= dict(
    include_vars=["longitude", "latitude"], indexs=dict(obs=slice(0,300))
)

filename="cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-0.25deg_P1D_1681265533678.nc"
eddies_collections = EddiesObservations.load_file(filename,**kwargs)

Actual outcome

 eddies_collections = EddiesObservations.load_file(filename,**kwargs)
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\observations\observation.py", line 782, in load_file
    return cls.load_from_netcdf(filename, **kwargs)
  File "C:\Users\dell\AppData\Local\Programs\Python\Python38\lib\site-packages\py_eddy_tracker\observations\observation.py", line 969, in load_from_netcdf
    nb_obs = len(h_nc.dimensions[obs_dim])
KeyError: None

It seems that this is not the problem of my file. I don't know why the error was reported. The following is the content of my file

<xarray.Dataset>
Dimensions:    (time: 5, latitude: 5, longitude: 5)
Coordinates:
  * latitude   (latitude) float32 3.875 4.125 4.375 4.625 4.875
  * time       (time) datetime64[ns] 2012-01-01 2012-01-02 ... 2012-01-05
  * longitude  (longitude) float32 104.9 105.1 105.4 105.6 105.9
Data variables:
    ugos       (time, latitude, longitude) float64 ...
    vgos       (time, latitude, longitude) float64 ...
    vgosa      (time, latitude, longitude) float64 ...
    crs        int32 ...
    err_vgosa  (time, latitude, longitude) float64 ...
    sla        (time, latitude, longitude) float64 ...
    ugosa      (time, latitude, longitude) float64 ...
    flag_ice   (time, latitude, longitude) float64 ...
    adt        (time, latitude, longitude) float64 ...
    err_ugosa  (time, latitude, longitude) float64 ...
    err_sla    (time, latitude, longitude) float64 ...
Attributes: (12/45)
    Conventions:                                    CF-1.6
    FROM_ORIGINAL_FILE__Metadata_Conventions:       Unidata Dataset Discovery...
    cdm_data_type:                                  Grid
    comment:                                        Sea Surface Height measur...
    contact:                                        servicedesk.cmems@mercato...
    creator_email:                                  servicedesk.cmems@mercato...
    ...                                             ...
    time_coverage_duration:                         P1D
    time_coverage_end:                              2022-06-23T12:00:00Z
    time_coverage_resolution:                       P1D
    time_coverage_start:                            2022-06-22T12:00:00Z
    title:                                          DT merged all satellites ...
    _CoordSysBuilder:                               ucar.nc2.dataset.conv.CF1...
AntSimi commented 1 year ago

It seem you try to load a gridded dataset with tools to load eddy dataset.