sciserver / poseidon-viewer

Apache License 2.0
3 stars 3 forks source link

ECCO interpolator error at level 4 #14

Open MaceKuailv opened 1 month ago

MaceKuailv commented 1 month ago

There is an issue when trying to find the 'G' node of a given lat lon, the new face and the original face is not connected, even when transitive is enabled. The errors are

--------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 1
----> 1 store_interpolator(
      2     env, zooms, subocedata, unique_grain, inverse_grain, exmp_vel=uv, exmp_scl=s
      3 )

File ~/workspace/Temporary/Thomas.Haine/scratch/poseidon-viewer/precalc/gen_interp.py:462, in store_interpolator(env, zooms, subocedata, unique_grains, inverse_grain, exmp_vel, exmp_scl)
    460 sub2use = subocedata[inverse_grain[iz]]
    461 grain = unique_grains[inverse_grain[iz]]
--> 462 interpolator = weight_index_inverse_from_latlon(
    463     sub2use,
    464     y,
    465     x,
    466     var=var,
    467     grain=grain,
    468     exmp_vel=exmp_vel,
    469     exmp_scl=exmp_scl,
    470 )
    471 txn.put(key, pickle.dumps(interpolator))

File ~/workspace/Temporary/Thomas.Haine/scratch/poseidon-viewer/precalc/gen_interp.py:416, in weight_index_inverse_from_latlon(oce, lat, lon, var, grain, exmp_vel, exmp_scl)
    382 def weight_index_inverse_from_latlon(
    383     oce, lat, lon, var="scalar", grain=None, exmp_vel=None, exmp_scl=None
    384 ):
    385     """Create the interpolator
    386 
    387     Parameters
   (...)
    414         If tuple, the second element (int) devides which part is for U and V,respectively.
    415     """
--> 416     pt = sd_position_from_latlon(lat, lon, oce)
    417     if var == "scalar":
    418         weight = calc_scalar_weight(pt)

File ~/workspace/Temporary/Thomas.Haine/scratch/poseidon-viewer/precalc/gen_interp.py:227, in sd_position_from_latlon(lat, lon, ocedata)
    225 need_rot = np.where(pt.face != pt.fcg)[0]
    226 for i in need_rot:
--> 227     edge, new_edge = tp.mutual_direction(pt.face[i], pt.fcg[i], transitive=True)
    228     rot = (
    229         np.pi - sd.topology.DIRECTIONS[edge] + sd.topology.DIRECTIONS[new_edge]
    230     ) % (np.pi * 2)
    231     pt.rxg[i], pt.ryg[i] = sd.utils.local_to_latlon(
    232         pt.rxg[i], pt.ryg[i], np.cos(rot), np.sin(rot)
    233     )

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/seaduck/topology.py:364, in Topology.mutual_direction(self, face, new_face, **kwarg)
    356 """Find the relative orientation of two faces.
    357 
    358 0,1,2,3 stands for up, down, left, right
   (...)
    361 2. the 1st face is to which direction of the 2nd face.
    362 """
    363 if self.typ == "LLC":
--> 364     return _llc_mutual_direction(face, new_face, **kwarg)
    365 elif self.typ in ["x_periodic", "box"]:
    366     raise ValueError(
    367         "It makes no sense to tinker with face_connection when there is only one face"
    368     )

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/seaduck/topology.py:55, in _llc_mutual_direction()
     53         break
     54 if common < 0:
---> 55     raise ValueError(
     56         "The two faces does not share common face, transitive did not help"
     57     )
     58 else:
     59     edge_1 = np.where(LLC_FACE_CONNECT[face] == common)[0][0]

ValueError: The two faces does not share common face, transitive did not help
MaceKuailv commented 1 month ago

This is a low priority issue because it must have resulted from the extreme mismatch between the grid spacing.