localdevices / pyorc

Surface velocity, object tracking, and river flow measurements in an open-source API
GNU Affero General Public License v3.0
129 stars 31 forks source link

Question about pyorc Tuturial | Transect of velocity values from the velocimetry results. #102

Closed moroth89 closed 1 year ago

moroth89 commented 1 year ago

Hi, I am doing the "quick start" from the pyorc 0.3.0 documentation (https://localdevices.github.io/pyorc/_examples/04_Extracting_crosssection_velocities_and_discharge.html). At the second last step "Obtain a discharge measurement over a cross section" I get an error as soon as I enter the fourth operation. I am using the latest pyopenrivercam version, and python 3.10.6.

In [4]: ds_points = ds.velocimetry.get_transect(x, y, z, crs=32735, rolling=4)^M
   ...: ds_points2 = ds.velocimetry.get_transect(x2, y2, z2, crs=32735, rolling=4)^M
   ...: ds_points
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 ds_points = ds.velocimetry.get_transect(x, y, z, crs=32735, rolling=4)
      2 ds_points2 = ds.velocimetry.get_transect(x2, y2, z2, crs=32735, rolling=4)
      3 ds_points

File ~\AppData\Local\mambaforge\lib\site-packages\pyorc\api\velocimetry.py:600, in Velocimetry.get_transect(self, x, y, z, s, crs, v_eff, xs, ys, distance, wdw, wdw_x_min, wdw_x_max, wdw_y_min, wdw_y_max, rolling, tolerance, quantiles)
    598     # remove velocities that are too few in samples
    599     ds_effective = ds_effective.where(missing_tolerance)
--> 600     ds_points = ds_effective.interp(x=_x, y=_y)
    601 if np.isnan(ds_points["v_x"].mean(dim="time")).all():
    602     warnings.warn(
    603         "No valid velocimetry points found over bathymetry. Check if the bathymetry is within the camera objective or anything is visible in objective."
    604     )

File ~\AppData\Local\mambaforge\lib\site-packages\xarray\core\dataset.py:3378, in Dataset.interp(self, coords, method, assume_sorted, kwargs, method_non_numeric, **coords_kwargs)
   3375 if dtype_kind in "uifc":
   3376     # For normal number types do the interpolation:
   3377     var_indexers = {k: v for k, v in use_indexers.items() if k in var.dims}
-> 3378     variables[name] = missing.interp(var, var_indexers, method, **kwargs)
   3379 elif dtype_kind in "ObU" and (use_indexers.keys() & var.dims):
   3380     # For types that we do not understand do stepwise
   3381     # interpolation to avoid modifying the elements.
   (...)
   3384     # this loop there might be some duplicate code that slows it
   3385     # down, therefore collect these signals and run it later:
   3386     reindex = True

File ~\AppData\Local\mambaforge\lib\site-packages\xarray\core\missing.py:639, in interp(var, indexes_coords, method, **kwargs)
    637 original_dims = broadcast_dims + dims
    638 new_dims = broadcast_dims + list(destination[0].dims)
--> 639 interped = interp_func(
    640     var.transpose(*original_dims).data, x, destination, method, kwargs
    641 )
    643 result = Variable(new_dims, interped, attrs=var.attrs)
    645 # dimension of the output array

File ~\AppData\Local\mambaforge\lib\site-packages\xarray\core\missing.py:764, in interp_func(var, x, new_x, method, kwargs)
    748         meta = var._meta
    750     return da.blockwise(
    751         _dask_aware_interpnd,
    752         out_ind,
   (...)
    761         align_arrays=False,
    762     )
--> 764 return _interpnd(var, x, new_x, func, kwargs)

File ~\AppData\Local\mambaforge\lib\site-packages\xarray\core\missing.py:788, in _interpnd(var, x, new_x, func, kwargs)
    786 # stack new_x to 1 vector, with reshape
    787 xi = np.stack([x1.values.ravel() for x1 in new_x], axis=-1)
--> 788 rslt = func(x, var, xi, **kwargs)
    789 # move back the interpolation axes to the last position
    790 rslt = rslt.transpose(range(-rslt.ndim + 1, 1))

File ~\AppData\Local\mambaforge\lib\site-packages\scipy\interpolate\_rgi.py:654, in interpn(points, values, xi, method, bounds_error, fill_value)
    650 if method in ["linear", "nearest", "slinear", "cubic", "quintic", "pchip"]:
    651     interp = RegularGridInterpolator(points, values, method=method,
    652                                      bounds_error=bounds_error,
    653                                      fill_value=fill_value)
--> 654     return interp(xi)
    655 elif method == "splinef2d":
    656     xi_shape = xi.shape

File ~\AppData\Local\mambaforge\lib\site-packages\scipy\interpolate\_rgi.py:336, in RegularGridInterpolator.__call__(self, xi, method)
    332 if (ndim == 2 and hasattr(self.values, 'dtype') and
    333         self.values.ndim == 2):
    334     # a fast path
    335     out = np.empty(indices.shape[1], dtype=self.values.dtype)
--> 336     result = evaluate_linear_2d(self.values,
    337                                 indices,
    338                                 norm_distances,
    339                                 self.grid,
    340                                 out)
    341 else:
    342     result = self._evaluate_linear(indices, norm_distances)

File _rgi_cython.pyx:19, in scipy.interpolate._rgi_cython.__pyx_fused_cpdef()

TypeError: No matching signature found

Up to this point, I was able to perform all steps as described in the "quick start". Also my cross sections are in the area of interest. Is it possible that there is a bug?

Thx a lot in advance for your help.

Cheers, Mo

hcwinsemius commented 1 year ago

Hi @moroth89 . Are you using our example or your own data? If you use your own could you send your input netcdf (I.e. results from the get_piv step)?

If not, could you send the output of pyorc.version? I can check this tomorrow.

hcwinsemius commented 1 year ago

Kindly also provide the x, y, z points of your transect and the coordinate reference system in which they have been measured if this is relevant.

moroth89 commented 1 year ago

Hi hcwinsemius,

I used the data from your example (video, coordinate reference system, gcps, etc.).
Here are the generated input ntcdf and the file from the get_piv step. ngwerere_mo.zip

Thanks a lot.

Cheers, Mo

hcwinsemius commented 1 year ago

Very useful that you noticed this, so thanks. The bug has been caused by a new release (yesterday!!) of scipy, causing problems in the interpolation routine of xarray. If have solved this for now in release v0.3.2. The expected rebuild on conda-forge is due tomorrow. Once it is there please update your package with:

conda update -c conda-forge pyopenrivercam

moroth89 commented 1 year ago

@hcwinsemius, Thanks a lot for your quick reply! It works, but I get different results. As you can see, some information is lost once I apply the filter in step 4.

Missing_values_along_cross-sections Missing_velocities_after_filtering

Please let me know if you would like me to open a new thread for this issue.

Thanks a lot in advance.

Cheers, Moritz

hcwinsemius commented 1 year ago

Moritz,

I am making significant changes in the filtering methods which cause that too many velocities are removed in your notebook. I am currently also describing the idea of the filters (which will be called "masks" in upcoming changes) and how to logically apply them. The notebook currently does not yet reflect that change.

I already have a changed recipe for that notebook which I can provide soon, probably next week when I have returned to office.

Yours sincerely.

Op ma 9 jan. 2023 13:03 schreef Moritz @.***>:

Thanks a lot for your quick reply! It works, but I get different results. As you can see, some information is lost once I apply the filter in point 4. https://localdevices.github.io/pyorc/_examples/03_Plotting_and_filtering_velocimetry_results.html [image: Missing_values_along_cross-sections] https://user-images.githubusercontent.com/59596749/211303021-7ebc6caa-0854-48a2-bd2b-64d7b331f945.PNG [image: Missing_velocities_after_filtering] https://user-images.githubusercontent.com/59596749/211303025-63a242f3-edab-40cb-a5a5-72c7d159908c.PNG

Please let me know if you would like me to open a new thread for this issue.

Thanks a lot in advance.

Cheers, Moritz

— Reply to this email directly, view it on GitHub https://github.com/localdevices/pyorc/issues/102#issuecomment-1375523302, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB2NZMP7B2S3SO2MUMC6PYLWRP5AXANCNFSM6AAAAAATQ7LYEY . You are receiving this because you modified the open/close state.Message ID: @.***>

moroth89 commented 1 year ago

Okay, I'm really curious about the new changes and look forward to digging deeper into this topic!

hcwinsemius commented 1 year ago

Just checked. If you want to use the very latest changes then you should check out branch "issue_101". I checked and it already contains the updated notebook (starting with 03...) with the explanation of the filters (aka masks). But you then also need to remove the conda installed package. If you wish to try that today then do the following within the folder co taining the git repository after activating your conda environment.

1 remove the conda-forge version

conda remove --force pyopenrivercam

2 checkout the right branch

git pull git checkout issue_101

3 install directly from code base

pip install -e .

Then try the new notebook with masks. I expect to release these changes somewhere next week. So you can also wait for that.

Op di 10 jan. 2023 09:20 schreef Moritz @.***>:

Okay, I'm really curious about the new changes and look forward to digging deeper into this topic!

— Reply to this email directly, view it on GitHub https://github.com/localdevices/pyorc/issues/102#issuecomment-1376888274, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB2NZMK7XF2ORJ23HKOSSQDWRULWFANCNFSM6AAAAAATQ7LYEY . You are receiving this because you were mentioned.Message ID: @.***>