Polygon-Based Data Cropping

[AtiehAlipour-NOAA]

Hi, I have a question: do we need to adjust this part of the code to use a polygon like the one below for cropping the data? Thanks!

polygon = array ([[271.09419606, 30.29241252],
[270.80483714, 30.47217758],
[270.31435613, 30.35380055],
[269.89991289, 30.0740687 ],
[269.9734827 , 29.72115943],
[270.38302128, 29.34776829],
[270.5129946 , 28.86997287],
[270.96177041, 28.87641531],
[271.38357063, 29.16376929],
[271.45468811, 29.68281716],
[271.68765915, 29.83609834],
[271.65087425, 29.96790232],
[271.09419606, 30.29241252]])

[pmav99]

You will need to update the implementation too. As it is it will only work with rectangular regions. For a polygon you will need to use some shapely method like https://shapely.readthedocs.io/en/stable/reference/shapely.contains.html or geopandas sjoin: https://geopandas.org/en/stable/docs/reference/api/geopandas.sjoin.html.

And you will need to check performance for large meshes. The pure numpy implementation we have here is probably faster than the alternatives. Nevertheless, try to get something working and if performance proves to be problematic you can investigate making it faster afterwards (e.g. leverage numba or something).

[josephzhang8]

I cannot find Searvey anymroe; anyone? Thx.

[AtiehAlipour-NOAA]

You will need to update the implementation too. As it is it will only work with rectangular regions. For a polygon you will need to use some shapely method like https://shapely.readthedocs.io/en/stable/reference/shapely.contains.html or geopandas sjoin: https://geopandas.org/en/stable/docs/reference/api/geopandas.sjoin.html.

And you will need to check performance for large meshes. The pure numpy implementation we have here is probably faster than the alternatives. Nevertheless, try to get something working and if performance proves to be problematic you can investigate making it faster afterwards (e.g. leverage numba or something).

Thanks for the quick response and the detailed explanation. For now, I am using a rectangle region but I wanted to confirm if there's a possibility to switch to a polygon later on. Thanks once again.

[AtiehAlipour-NOAA]

I cannot find Searvey anymroe; anyone? Thx.

Hi, please check here: https://github.com/oceanmodeling/searvey

[josephzhang8]

Thx, Atieh!

[pmav99]

https://github.com/asascience-open/xarray-subset-grid

[AtiehAlipour-NOAA]

I updated the crop function in Thalassa/utils.py to subset stofs-3d data using a polygon:

def crop(
    ds: xarray.Dataset,
    bbox: shapely.Polygon,
) -> xarray.Dataset:
    """
    Crop the dataset using the provided `bbox`, which can be a bounding box or an arbitrary polygon.

    Examples:
        ``` python
        import thalassa
        import shapely
        from shapely.geometry import Polygon
        
        ds = thalassa.open_dataset("some_netcdf.nc")
        bbox = shapely.box(0, 0, 1, 1) 
        or 
        bbox = Polygon([
        (10, 40),
        (20, 40),
        (15, 50)
        ])
        ds = thalassa.crop(ds, bbox)
        ```

    Parameters:
        ds: The dataset we want to crop.
        bbox: A Shapely polygon whose boundary will be used to crop `ds`.
    """

 
    import numpy as np
    import numpy_indexed as npi
    import shapely.geometry

    bbox = resolve_bbox(bbox)

    node_points = [shapely.geometry.Point(lon, lat) for lon, lat in zip(ds.lon.values, ds.lat.values)]

    mask = np.array([bbox.contains(point) for point in node_points])

    indices_of_nodes_in_bbox = np.where(mask)[0]

    triface_nodes = ds.triface_nodes.values
    triangle_mask = np.all(np.isin(triface_nodes, indices_of_nodes_in_bbox), axis=1)
    indices_of_triface_nodes_in_bbox = np.where(triangle_mask)[0]

    ds = ds.isel(node=indices_of_nodes_in_bbox, triface=indices_of_triface_nodes_in_bbox)

    remapped_nodes = np.arange(len(indices_of_nodes_in_bbox))
    remapped_triface_nodes = np.c_[
        npi.remap(ds.triface_nodes[:, 0], indices_of_nodes_in_bbox, remapped_nodes),
        npi.remap(ds.triface_nodes[:, 1], indices_of_nodes_in_bbox, remapped_nodes),
        npi.remap(ds.triface_nodes[:, 2], indices_of_nodes_in_bbox, remapped_nodes),
    ]
    ds["triface_nodes"] = (("triface", "three"), remapped_triface_nodes)

    return ds

and it worked: