tobac-project / tobac

Tracking and object-based analysis of clouds
BSD 3-Clause "New" or "Revised" License
99 stars 53 forks source link

Problem with monotonic longitude requirements #382

Open wreckdump opened 9 months ago

wreckdump commented 9 months ago

I would like to use tobac with the input that has non-monotonic longitude. For example, below attached picture shows the entire tropics (25S - 25N) with the longitude ranging from 30E - 27.75E, which means that in the middle of the Pacific Ocean the longitude's coordinate transitions from 180E to 180W, and so on. This makes the longitude non-monotonic.

scrshot 2023 12 06 14 40 45

How can I work around this issue?

wreckdump commented 9 months ago

Aha! Thank you so much! Just so happens to be that I am using version 1.5.1! I'll try this! Thank you again!

But, would I be still getting the error message of

/home/rangke/py.vrt/lib/python3.11/site-packages/iris/fileformats/_nc_load_rules/helpers.py:1085: UserWarning: Failed to create 'lon' dimension coordinate: The 'longitude' DimCoord points array must be strictly monotonic.
Gracefully creating 'lon' auxiliary coordinate instead.
  warnings.warn(msg.format(name=str(cf_coord_var.cf_name), error=e_msg))

even when I specify PCB_flag in the tobac.feature_detection.feature_detection_multithreshold()?

Thanks.

freemansw1 commented 9 months ago

Deleted my original comment as I misunderstood the question. The PBC_flag will still be necessary here as you have wrap-around data though.

In theory, this should work fine without doing anything, but I'm guessing there could be an issue with add_coords. The simple answer would be to interpolate the coordinates yourself, but we can flag this and take a look.

For your specific question, that warning appears to be from iris, and can perhaps be ignored? As I mentioned, you may have an issue with the coordinates you get out of tobac feature detection, though. Take a look at the longitude output and see if it gives odd answers near the +/- 180 boundary

w-k-jones commented 9 months ago

To add to what @freemansw1 says, most of tobac will run fine with non-monotonic coordinates, however area calculation will fail and the latitude/longitude of features that cross the date line will be wrong due to using linear interpolation (tracking/segmentation will still work correctly because these use array index rather than lat/lon).

An easy fix to these issues would be to load your data in xarray and roll the data so that the longitude is monotonic, then convert to iris for use with tobac. e.g.:

import xarray as xr
bt = xr.DataArray("path/to/data.nc")
bt = bt.roll({"lon":offset}, roll_coords=True)
bt = bt.load().to_iris()

features = tobac.feature_detection(bt, ...)

The required offset is in the index values, so if your data has 0.1 degree grid spacing you would need to roll by 2100 to get the lon array to span -180->180.

Note also that PBCs don't work with predictive tracking in v1.5.1, but this is fixed for v1.5.2 which will be released very soon

wreckdump commented 9 months ago

Errm. I don't think the periodic boundary condition function is working properly. I was just testing out a sample data that I have made (https://file.io/XjmoIGLdXurn) with the following info.

dimensions:
    time = UNLIMITED ; // (22 currently)
    lat = 201 ;
    lon = 1440 ;
variables:
    float time(time) ;
        time:standard_name = "time" ;
        time:units = "hours since 2020-01-21 00:00:00" ;
    float lat(lat) ;
        lat:standard_name = "latitude" ;
        lat:units = "degree_north" ;
    float lon(lon) ;
        lon:standard_name = "longitude" ;
        lon:units = "degree_east" ;
    float tb(time, lat, lon) ;
        tb:_FillValue = 99999.f ;
        tb:standard_name = "brightness temperature" ;
        tb:units = "K" ;

And some of the snapshots of the file looks like below. Longitude is from 180W to 180E (or -180 to 180), which is monotonic.

scrshot 2023 12 06 17 21 17 scrshot 2023 12 06 17 21 35 scrshot 2023 12 06 17 21 51

Background value is 300, and those square patches are set to 230. Tracked results are shown below.

trk pnc err

I expected those two different colored tracks to be the same trajectory group but they are separated in two. Here are the options that I have used.

   # feature detection
   param_fet={}
   param_fet['target']='minimum'
   param_fet['threshold']=[thres]
   param_fet['n_min_threshold']=4
   param_fet['position_threshold']='center'
   param_fet['sigma_threshold']=0.5
   param_fet['PBC_flag'] = 'hdim_2'

   # segmentation
   param_seg={}
   param_seg['target']='minimum'
   param_seg['method']='watershed'
   param_seg['threshold']=thres

  # linking 
   param_lnk={}
   param_lnk['v_max']=None
   param_lnk['d_max']=111.0e3*1 # 10 degrees in meters
   param_lnk['stubs']=1
   param_lnk['order']=2
   param_lnk['extrapolate']=0
   param_lnk['memory']=0
   param_lnk['adaptive_stop']=0.5
   param_lnk['adaptive_step']=0.95
   param_lnk['subnetwork_size']=10
   param_lnk['method_linking']= 'random'

I tried both random and predict for method_linking.

w-k-jones commented 9 months ago

Errm. I don't think the periodic boundary condition function is working properly. I was just testing out a sample data that I have made (https://file.io/XjmoIGLdXurn) with the following info.

I think you can fix this by adding the PBC_flag keyword to segmentation and tracking, and also addting min_h2 and max_h2 for tracking as below:

   # feature detection
   param_fet={}
   param_fet['target']='minimum'
   param_fet['threshold']=[thres]
   param_fet['n_min_threshold']=4
   param_fet['position_threshold']='center'
   param_fet['sigma_threshold']=0.5
   param_fet['PBC_flag'] = 'hdim_2'

   # segmentation
   param_seg={}
   param_seg['target']='minimum'
   param_seg['method']='watershed'
   param_seg['threshold']=thres
   param_seg['PBC_flag'] = 'hdim_2'

  # linking 
   param_lnk={}
   param_lnk['v_max']=None
   param_lnk['d_max']=111.0e3*1 # 10 degrees in meters
   param_lnk['stubs']=1
   param_lnk['order']=2
   param_lnk['extrapolate']=0
   param_lnk['memory']=0
   param_lnk['adaptive_stop']=0.5
   param_lnk['adaptive_step']=0.95
   param_lnk['subnetwork_size']=10
   param_lnk['method_linking']= 'random'
   param_lnk['PBC_flag'] = 'hdim_2'
   param_lnk['min_h2'] = 0
   param_lnk['max_h2'] = 1440

Note that min_h2 and max_h2 specify the size of your input array in the PBC dimension

I tried both random and predict for method_linking.

If you want to try PBCs with predictive tracking right away, you could checkout the RC_v1.5.x branch of the repository, although it should be released on conda soon

This is useful insight on using PBCs, it would be nice to have a way of storing the PBC settings so that you don't have to repeat it...

wreckdump commented 9 months ago

Thank you for your reply. After making changes as you suggested, I am getting the following error messages.

Traceback (most recent call last):
  File "/home/rangke/work/thesis/py.scr/trk.fet.seg.lnk.py", line 128, in <module>
    trk               =  tobac.tracking.linking_trackpy(fet_var, tb, dt=dt, dxy=dxy, **param_lnk)
                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/tobac/tracking.py", line 343, in linking_trackpy
    trajectories_unfiltered = tp.link(
                              ^^^^^^^^
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/trackpy/linking/linking.py", line 192, in link
    for i, _ids in link_iter(coords_iter, search_range, **kwargs):
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/trackpy/linking/linking.py", line 98, in link_iter
    linker.init_level(coords, t)
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/trackpy/linking/linking.py", line 476, in init_level
    self.update_hash(coords, t, extra_data)
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/trackpy/linking/linking.py", line 461, in update_hash
    self.hash = self.hash_cls(points_from_arr(coords, t, extra_data),
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/trackpy/linking/subnet.py", line 167, in __init__
    self.rebuild()
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/trackpy/linking/subnet.py", line 199, in rebuild
    self._btree = BallTree(coords_mapped,
                  ^^^^^^^^^^^^^^^^^^^^^^^
  File "sklearn/neighbors/_binary_tree.pxi", line 860, in sklearn.neighbors._ball_tree.BinaryTree.__init__
  File "sklearn/neighbors/_binary_tree.pxi", line 1039, in sklearn.neighbors._ball_tree.BinaryTree._recursive_build
  File "sklearn/neighbors/_ball_tree.pyx", line 101, in sklearn.neighbors._ball_tree.init_node
  File "sklearn/neighbors/_binary_tree.pxi", line 1017, in sklearn.neighbors._ball_tree.BinaryTree.rdist
  File "sklearn/metrics/_dist_metrics.pyx", line 353, in sklearn.metrics._dist_metrics.DistanceMetric64.rdist
  File "sklearn/metrics/_dist_metrics.pyx", line 2640, in sklearn.metrics._dist_metrics.PyFuncDistance64.dist
  File "sklearn/metrics/_dist_metrics.pyx", line 2651, in sklearn.metrics._dist_metrics.PyFuncDistance64._dist
  File "/home/rangke/py.vrt/lib/python3.11/site-packages/numba/core/dispatcher.py", line 468, in _compile_for_args
    error_rewrite(e, 'typing')
  File "/home/rangke/py.vrt/lib/python3.11/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)
No implementation of function Function(<built-in function sub>) found for signature:

 >>> sub(none, none)

There are 14 candidate implementations:
  - Of which 10 did not match due to:
  Overload of function 'sub': File: <numerous>: Line N/A.
    With argument(s): '(none, none)':
   No match.
  - Of which 2 did not match due to:
  Operator Overload in function 'sub': File: unknown: Line unknown.
    With argument(s): '(none, none)':
   No match for registered cases:
    * (int64, int64) -> int64
    * (int64, uint64) -> int64
    * (uint64, int64) -> int64
    * (uint64, uint64) -> uint64
    * (float32, float32) -> float32
    * (float64, float64) -> float64
    * (complex64, complex64) -> complex64
    * (complex128, complex128) -> complex128
  - Of which 2 did not match due to:
  Overload in function 'impl_set_difference': File: numba/cpython/setobj.py: Line 1525.
    With argument(s): '(none, none)':
   Rejected as the implementation raised a specific error:
     TypingError: All arguments must be Sets, got (none, none)
  raised from /home/rangke/py.vrt/lib/python3.11/site-packages/numba/cpython/setobj.py:108

During: typing of intrinsic-call at /home/rangke/py.vrt/lib/python3.11/site-packages/tobac/utils/periodic_boundaries.py (243)

File "../../../py.vrt/lib/python3.11/site-packages/tobac/utils/periodic_boundaries.py", line 243:
def calc_distance_coords_pbc(
    <source elided>
    if PBC_flag in ["hdim_1", "both"]:
        size_h1 = max_h1 - min_h1
        ^

I am not sure what the problem is here....

w-k-jones commented 9 months ago

Are you running this in a jupyter notebook? I've encountered a similar bug this evening that running PBCs with numba which occurs in notebooks but not in the python interpreter

wreckdump commented 9 months ago

No. I am not using jupyter notebook. I am just using python 3.11.6 with all the modules installed via pip.

w-k-jones commented 9 months ago

Hmmm, interesting. I can recreate the same error in a notebook but it runs fine from the command line... I will try and fix this before we release v1.5.2. Could you let me know what version of numba you have installed?

In the mean time, you can get tracking to work by uninstalling numba. It will run slower, but shouldn't cause an error

wreckdump commented 9 months ago

Numba version that I am using is 0.58.1, and I am also running tobac from terminal without any development environment.

The code runs without numba fine. Thanks.

freemansw1 commented 9 months ago

This is useful insight on using PBCs, it would be nice to have a way of storing the PBC settings so that you don't have to repeat it...

I have thought this for a while. We could probably implement such a change with the xarray transition through the attributes.