From b774e6380c781703d277e397b0d51be6c44c7daf Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Wed, 23 Nov 2022 23:52:33 +0100 Subject: [PATCH 01/11] WIP: add xvec accessor and coords_to_crs --- xvec/__init__.py | 1 + xvec/accessor.py | 59 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+) create mode 100644 xvec/accessor.py diff --git a/xvec/__init__.py b/xvec/__init__.py index 15ad2cb..44c022b 100644 --- a/xvec/__init__.py +++ b/xvec/__init__.py @@ -1,5 +1,6 @@ from importlib.metadata import PackageNotFoundError, version +from .accessor import XvecAccessor # noqa from .index import GeometryIndex # noqa try: diff --git a/xvec/accessor.py b/xvec/accessor.py new file mode 100644 index 0000000..996eb7a --- /dev/null +++ b/xvec/accessor.py @@ -0,0 +1,59 @@ +from typing import Any, Hashable, Union + +import numpy as np +import shapely +import xarray as xr +from pyproj import CRS, Transformer + +from .index import GeometryIndex + + +@xr.register_dataarray_accessor("xvec") +@xr.register_dataset_accessor("xvec") +class XvecAccessor: + def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): + self._obj = xarray_obj + + def coords_to_crs(self, coords: Hashable, crs: Any): + + data = self._obj[coords] + data_crs = self._obj.xindexes[coords].crs + + # transformation code taken from geopandas (BSD 3-clause license) + if data_crs is None: + raise ValueError( + "Cannot transform naive geometries. " + "Please set a CRS on the object first." + ) + + crs = CRS.from_user_input(crs) + + if data_crs.is_exact_same(crs): + return self + + transformer = Transformer.from_crs(data_crs, crs, always_xy=True) + + has_z = shapely.has_z(data) + + result = np.empty_like(data) + + coordinates = shapely.get_coordinates(data[~has_z], include_z=False) + new_coords_z = transformer.transform(coordinates[:, 0], coordinates[:, 1]) + result[~has_z] = shapely.set_coordinates( + data[~has_z].copy(), np.array(new_coords_z).T + ) + + coords_z = shapely.get_coordinates(data[has_z], include_z=True) + new_coords_z = transformer.transform( + coords_z[:, 0], coords_z[:, 1], coords_z[:, 2] + ) + result[has_z] = shapely.set_coordinates( + data[has_z].copy(), np.array(new_coords_z).T + ) + + # return result + return ( + self._obj.assign_coords({coords: result}) + .drop_indexes(coords) + .set_xindex(coords, GeometryIndex, crs=crs) + ) From 394bc97ad443fb7009c5eec33c41637ae8115cb6 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Thu, 24 Nov 2022 00:00:33 +0100 Subject: [PATCH 02/11] fix the shortcut if crs is equal --- xvec/accessor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xvec/accessor.py b/xvec/accessor.py index 996eb7a..f092925 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -29,7 +29,7 @@ def coords_to_crs(self, coords: Hashable, crs: Any): crs = CRS.from_user_input(crs) if data_crs.is_exact_same(crs): - return self + return self._obj transformer = Transformer.from_crs(data_crs, crs, always_xy=True) From 88c837989db3560e0c39ebdcb35d82fea7277f4d Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Thu, 1 Dec 2022 23:10:26 +0100 Subject: [PATCH 03/11] adapt for standard syntax --- xvec/accessor.py | 90 ++++++++++++++++++++++++++++++------------------ 1 file changed, 56 insertions(+), 34 deletions(-) diff --git a/xvec/accessor.py b/xvec/accessor.py index f092925..a3be64b 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -1,4 +1,4 @@ -from typing import Any, Hashable, Union +from typing import Any, Mapping, Union import numpy as np import shapely @@ -14,46 +14,68 @@ class XvecAccessor: def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): self._obj = xarray_obj - def coords_to_crs(self, coords: Hashable, crs: Any): - - data = self._obj[coords] - data_crs = self._obj.xindexes[coords].crs - - # transformation code taken from geopandas (BSD 3-clause license) - if data_crs is None: + def coords_to_crs( + self, + coords: Mapping[Any, Any] | None = None, + **coords_kwargs: Any, + ): + if coords and coords_kwargs: raise ValueError( - "Cannot transform naive geometries. " - "Please set a CRS on the object first." + "cannot specify both keyword and positional arguments to " + ".xvec.coords_to_crs" ) - crs = CRS.from_user_input(crs) + _obj = self._obj.copy(deep=False) + + if coords_kwargs: + coords = coords_kwargs + + transformed = {} + + for key, crs in coords.items(): + + data = _obj[key] + data_crs = self._obj.xindexes[key].crs + + # transformation code taken from geopandas (BSD 3-clause license) + if data_crs is None: + raise ValueError( + "Cannot transform naive geometries. " + f"Please set a CRS on the '{key}' coordinates first." + ) - if data_crs.is_exact_same(crs): - return self._obj + crs = CRS.from_user_input(crs) - transformer = Transformer.from_crs(data_crs, crs, always_xy=True) + if data_crs.is_exact_same(crs): + pass - has_z = shapely.has_z(data) + transformer = Transformer.from_crs(data_crs, crs, always_xy=True) - result = np.empty_like(data) + has_z = shapely.has_z(data) - coordinates = shapely.get_coordinates(data[~has_z], include_z=False) - new_coords_z = transformer.transform(coordinates[:, 0], coordinates[:, 1]) - result[~has_z] = shapely.set_coordinates( - data[~has_z].copy(), np.array(new_coords_z).T - ) + result = np.empty_like(data) - coords_z = shapely.get_coordinates(data[has_z], include_z=True) - new_coords_z = transformer.transform( - coords_z[:, 0], coords_z[:, 1], coords_z[:, 2] - ) - result[has_z] = shapely.set_coordinates( - data[has_z].copy(), np.array(new_coords_z).T - ) + coordinates = shapely.get_coordinates(data[~has_z], include_z=False) + new_coords_z = transformer.transform(coordinates[:, 0], coordinates[:, 1]) + result[~has_z] = shapely.set_coordinates( + data[~has_z].copy(), np.array(new_coords_z).T + ) + + coords_z = shapely.get_coordinates(data[has_z], include_z=True) + new_coords_z = transformer.transform( + coords_z[:, 0], coords_z[:, 1], coords_z[:, 2] + ) + result[has_z] = shapely.set_coordinates( + data[has_z].copy(), np.array(new_coords_z).T + ) + + transformed[key] = (result, crs) + + for key, (result, crs) in transformed.items(): + _obj = ( + _obj.assign_coords({key: result}) + .drop_indexes(key) + .set_xindex(key, GeometryIndex, crs=crs) + ) - # return result - return ( - self._obj.assign_coords({coords: result}) - .drop_indexes(coords) - .set_xindex(coords, GeometryIndex, crs=crs) - ) + return _obj From 468fba88f669e12407c3899f2aea1087fce9ad5f Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Fri, 2 Dec 2022 14:41:59 +0100 Subject: [PATCH 04/11] implement the agreed API --- xvec/accessor.py | 92 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 77 insertions(+), 15 deletions(-) diff --git a/xvec/accessor.py b/xvec/accessor.py index a3be64b..1ee19f7 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -1,4 +1,4 @@ -from typing import Any, Mapping, Union +from typing import Any, List, Mapping, Union import numpy as np import shapely @@ -13,26 +13,35 @@ class XvecAccessor: def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): self._obj = xarray_obj + self._geom_coord_names = [ + name + for name, index in self._obj.xindexes.items() + if isinstance(index, GeometryIndex) + ] - def coords_to_crs( + @property + def geom_coords_names(self) -> List: + return self._geom_coord_names + + def to_crs( self, - coords: Mapping[Any, Any] | None = None, - **coords_kwargs: Any, + variable_crs: Mapping[Any, Any] | None = None, + **variable_crs_kwargs: Any, ): - if coords and coords_kwargs: + if variable_crs and variable_crs_kwargs: raise ValueError( "cannot specify both keyword and positional arguments to " - ".xvec.coords_to_crs" + ".xvec.to_crs" ) _obj = self._obj.copy(deep=False) - if coords_kwargs: - coords = coords_kwargs + if variable_crs_kwargs: + variable_crs = variable_crs_kwargs transformed = {} - for key, crs in coords.items(): + for key, crs in variable_crs.items(): data = _obj[key] data_crs = self._obj.xindexes[key].crs @@ -56,9 +65,9 @@ def coords_to_crs( result = np.empty_like(data) coordinates = shapely.get_coordinates(data[~has_z], include_z=False) - new_coords_z = transformer.transform(coordinates[:, 0], coordinates[:, 1]) + new_coords = transformer.transform(coordinates[:, 0], coordinates[:, 1]) result[~has_z] = shapely.set_coordinates( - data[~has_z].copy(), np.array(new_coords_z).T + data[~has_z].copy(), np.array(new_coords).T ) coords_z = shapely.get_coordinates(data[has_z], include_z=True) @@ -72,10 +81,63 @@ def coords_to_crs( transformed[key] = (result, crs) for key, (result, crs) in transformed.items(): - _obj = ( - _obj.assign_coords({key: result}) - .drop_indexes(key) - .set_xindex(key, GeometryIndex, crs=crs) + _obj = _obj.assign_coords({key: result}) + + # TODO: use this instead of what is below once pydata/xarray#7347 is released + # _obj = _obj.drop_indexes(variable_crs.keys()) + + # for key, crs in variable_crs.items(): + # _obj = _obj.set_xindex(key, GeometryIndex, crs=crs) + + # this until return is a workaround for pydata/xarray#7347 + unchanged_geom_coords = { + name: self._obj.xindexes[name].crs + for name in self.geom_coords_names + if name not in variable_crs + } + + all_geom_coords = {**unchanged_geom_coords, **variable_crs} + + _obj = _obj.drop_indexes(all_geom_coords.keys()) + + for key, crs in all_geom_coords.items(): + _obj = _obj.set_xindex(key, GeometryIndex, crs=crs) + + return _obj + + def set_crs( + self, + variable_crs: Mapping[Any, Any] | None = None, + allow_override=False, + **variable_crs_kwargs: Any, + ): + + if variable_crs and variable_crs_kwargs: + raise ValueError( + "cannot specify both keyword and positional arguments to " + ".xvec.set_crs" ) + _obj = self._obj.copy(deep=False) + + if variable_crs_kwargs: + variable_crs = variable_crs_kwargs + + for key, crs in variable_crs.items(): + + data_crs = self._obj.xindexes[key].crs + + if not allow_override and data_crs is not None and not data_crs == crs: + raise ValueError( + f"The index '{key}' already has a CRS which is not equal to the " + "passed CRS. Specify 'allow_override=True' to allow replacing the " + "existing CRS without doing any transformation. If you actually " + "want to transform the geometries, use '.xvec.to_crs' instead." + ) + + _obj = _obj.drop_indexes(variable_crs_kwargs.keys()) + + for key, crs in variable_crs_kwargs.items(): + _obj = _obj.set_xindex(key, GeometryIndex, crs=crs) + return _obj From 6d6e1439c7b7ccb0839d3b9de0da5a4d891e31b3 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Sat, 3 Dec 2022 15:17:09 +0100 Subject: [PATCH 05/11] use the fix from xarray 2022.12 --- pyproject.toml | 2 +- xvec/accessor.py | 19 ++----------------- 2 files changed, 3 insertions(+), 18 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 1333517..c97c6fb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,7 +27,7 @@ classifiers = [ ] requires-python = ">=3.8" dependencies = [ - "xarray >= 2022.11.0", + "xarray >= 2022.12.0", "pyproj >= 3.0.0", "shapely >= 2.0b1", ] diff --git a/xvec/accessor.py b/xvec/accessor.py index 1ee19f7..7b51c13 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -83,24 +83,9 @@ def to_crs( for key, (result, crs) in transformed.items(): _obj = _obj.assign_coords({key: result}) - # TODO: use this instead of what is below once pydata/xarray#7347 is released - # _obj = _obj.drop_indexes(variable_crs.keys()) + _obj = _obj.drop_indexes(variable_crs.keys()) - # for key, crs in variable_crs.items(): - # _obj = _obj.set_xindex(key, GeometryIndex, crs=crs) - - # this until return is a workaround for pydata/xarray#7347 - unchanged_geom_coords = { - name: self._obj.xindexes[name].crs - for name in self.geom_coords_names - if name not in variable_crs - } - - all_geom_coords = {**unchanged_geom_coords, **variable_crs} - - _obj = _obj.drop_indexes(all_geom_coords.keys()) - - for key, crs in all_geom_coords.items(): + for key, crs in variable_crs.items(): _obj = _obj.set_xindex(key, GeometryIndex, crs=crs) return _obj From 543e76a1f824497a3525a913314b9685288bd266 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Sat, 3 Dec 2022 16:54:55 +0100 Subject: [PATCH 06/11] test coverage --- xvec/accessor.py | 16 ++-- xvec/tests/conftest.py | 61 +++++++++++++++ xvec/tests/test_accessor.py | 146 ++++++++++++++++++++++++++++++++++++ xvec/tests/test_index.py | 34 --------- 4 files changed, 215 insertions(+), 42 deletions(-) create mode 100644 xvec/tests/conftest.py create mode 100644 xvec/tests/test_accessor.py diff --git a/xvec/accessor.py b/xvec/accessor.py index 7b51c13..3e09c08 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -13,7 +13,7 @@ class XvecAccessor: def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): self._obj = xarray_obj - self._geom_coord_names = [ + self._geom_coords_names = [ name for name, index in self._obj.xindexes.items() if isinstance(index, GeometryIndex) @@ -21,7 +21,7 @@ def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): @property def geom_coords_names(self) -> List: - return self._geom_coord_names + return self._geom_coords_names def to_crs( self, @@ -30,8 +30,8 @@ def to_crs( ): if variable_crs and variable_crs_kwargs: raise ValueError( - "cannot specify both keyword and positional arguments to " - ".xvec.to_crs" + "Cannot specify both keyword and positional arguments to " + "'.xvec.to_crs'." ) _obj = self._obj.copy(deep=False) @@ -99,8 +99,8 @@ def set_crs( if variable_crs and variable_crs_kwargs: raise ValueError( - "cannot specify both keyword and positional arguments to " - ".xvec.set_crs" + "Cannot specify both keyword and positional arguments to " + ".xvec.set_crs." ) _obj = self._obj.copy(deep=False) @@ -120,9 +120,9 @@ def set_crs( "want to transform the geometries, use '.xvec.to_crs' instead." ) - _obj = _obj.drop_indexes(variable_crs_kwargs.keys()) + _obj = _obj.drop_indexes(variable_crs.keys()) - for key, crs in variable_crs_kwargs.items(): + for key, crs in variable_crs.items(): _obj = _obj.set_xindex(key, GeometryIndex, crs=crs) return _obj diff --git a/xvec/tests/conftest.py b/xvec/tests/conftest.py new file mode 100644 index 0000000..c67b9e8 --- /dev/null +++ b/xvec/tests/conftest.py @@ -0,0 +1,61 @@ +import numpy as np +import pytest +import shapely +import xarray as xr +from pyproj import CRS + +from xvec import GeometryIndex + + +@pytest.fixture(scope="session") +def geom_array(): + return np.array([shapely.Point(1, 2), shapely.Point(3, 4)]) + + +@pytest.fixture(scope="session") +def geom_array_z(): + return np.array([shapely.Point(10, 20, 30), shapely.Point(30, 40, 50)]) + + +@pytest.fixture(scope="session") +def geom_dataset_no_index(geom_array): + # a dataset with a geometry coordinate but no index + ds = xr.Dataset(coords={"geom": geom_array}) + return ds.drop_indexes("geom") + + +@pytest.fixture(scope="session") +def geom_dataset(geom_dataset_no_index): + # a dataset with a geometry coordinate baked by a GeometryIndex + crs = CRS.from_user_input(26915) + return geom_dataset_no_index.set_xindex("geom", GeometryIndex, crs=crs) + + +@pytest.fixture(scope="session") +def geom_dataset_no_crs(geom_dataset_no_index): + # a dataset with a geometry coordinate baked by a GeometryIndex (no CRS) + return geom_dataset_no_index.set_xindex("geom", GeometryIndex) + + +@pytest.fixture(scope="session") +def first_geom_dataset(geom_dataset, geom_array): + return ( + xr.Dataset(coords={"geom": [geom_array[0]]}) + .drop_indexes("geom") + .set_xindex("geom", GeometryIndex, crs=geom_dataset.xindexes["geom"].crs) + ) + + +@pytest.fixture(scope="session") +def multi_geom_dataset(geom_array, geom_array_z): + return ( + xr.Dataset( + coords={ + "geom": geom_array, + "geom_z": geom_array_z, + } + ) + .drop_indexes(["geom", "geom_z"]) + .set_xindex("geom", GeometryIndex, crs=26915) + .set_xindex("geom_z", GeometryIndex, crs=26915) + ) diff --git a/xvec/tests/test_accessor.py b/xvec/tests/test_accessor.py new file mode 100644 index 0000000..76d95a2 --- /dev/null +++ b/xvec/tests/test_accessor.py @@ -0,0 +1,146 @@ +import numpy as np +import pytest +import shapely +import xarray as xr + +import xvec # noqa +from xvec import GeometryIndex + + +@pytest.fixture() +def geom_array_4326(): + return shapely.from_wkt( + ["POINT (-97.488735 0.000018)", "POINT (-97.488717 0.000036)"] + ) + + +@pytest.fixture() +def geom_array_z_4326(): + return shapely.from_wkt( + ["POINT Z (-97.488654 0.00018 30)", "POINT Z (-97.488475 0.000361 50)"] + ) + + +@pytest.fixture() +def mixed_geom_dims_dataset(geom_array, geom_array_z): + return ( + xr.Dataset(coords={"geom": np.concatenate([geom_array, geom_array_z])}) + .drop_indexes("geom") + .set_xindex("geom", GeometryIndex, crs=26915) + ) + + +@pytest.fixture() +def geom_array_mixed_4326(geom_array_4326, geom_array_z_4326): + return np.concatenate([geom_array_4326, geom_array_z_4326]) + + +# Test .xvec accessor + + +def test_accessor(multi_geom_dataset): + assert hasattr(multi_geom_dataset, "xvec") + xr.testing.assert_identical(multi_geom_dataset.xvec._obj, multi_geom_dataset) + + +# Test .xvec geom_coords_names + + +def test_geom_coords_names(multi_geom_dataset): + assert multi_geom_dataset.xvec._geom_coords_names == ["geom", "geom_z"] + assert multi_geom_dataset.xvec.geom_coords_names == ["geom", "geom_z"] + + with pytest.raises(AttributeError, match="can't set attribute 'geom_coords_names'"): + multi_geom_dataset.xvec.geom_coords_names = ["a", "b"] + + +# Test .xvec.to_crs + + +def test_to_crs(geom_dataset, geom_array_4326): + transformed = geom_dataset.xvec.to_crs(geom=4326) + assert shapely.equals_exact(geom_array_4326, transformed.geom, tolerance=1e-5).all() + assert transformed.xindexes["geom"].crs.equals(4326) + + +def test_to_crs_two_geoms(multi_geom_dataset, geom_array_4326, geom_array_z_4326): + transformed = multi_geom_dataset.xvec.to_crs(geom=4326, geom_z=4326) + assert shapely.equals_exact(geom_array_4326, transformed.geom, tolerance=1e-5).all() + assert shapely.equals_exact( + geom_array_z_4326, transformed.geom_z, tolerance=1e-5 + ).all() + assert transformed.xindexes["geom"].crs.equals(4326) + assert transformed.xindexes["geom_z"].crs.equals(4326) + + +def test_to_crs_dict(multi_geom_dataset, geom_array_4326, geom_array_z_4326): + transformed = multi_geom_dataset.xvec.to_crs({"geom": 4326, "geom_z": 4326}) + assert shapely.equals_exact(geom_array_4326, transformed.geom, tolerance=1e-5).all() + assert shapely.equals_exact( + geom_array_z_4326, transformed.geom_z, tolerance=1e-5 + ).all() + assert transformed.xindexes["geom"].crs.equals(4326) + assert transformed.xindexes["geom_z"].crs.equals(4326) + + +def test_to_crs_arg_err(multi_geom_dataset): + with pytest.raises(ValueError, match="Cannot specify both keyword"): + multi_geom_dataset.xvec.to_crs(geom=4326, variable_crs={"geom_z": 4326}) + + +def test_to_crs_naive(geom_dataset_no_crs): + with pytest.raises(ValueError, match="Cannot transform naive geometries."): + geom_dataset_no_crs.xvec.to_crs(geom=4326) + + +def test_to_crs_same(geom_dataset): + transformed = geom_dataset.xvec.to_crs(geom=26915) + xr.testing.assert_identical(transformed, geom_dataset) + + +def test_to_crs_mixed_geom_dims(mixed_geom_dims_dataset, geom_array_mixed_4326): + transformed = mixed_geom_dims_dataset.xvec.to_crs(geom=4326) + assert shapely.equals_exact( + geom_array_mixed_4326, transformed.geom, tolerance=1e-5 + ).all() + assert transformed.xindexes["geom"].crs.equals(4326) + + +# Test .xvec.set_crs + + +def test_set_crs(geom_dataset_no_crs): + with_crs = geom_dataset_no_crs.xvec.set_crs(geom=4326) + assert with_crs.xindexes["geom"].crs.equals(4326) + + +def test_set_crs_dict(geom_dataset_no_crs): + with_crs = geom_dataset_no_crs.xvec.set_crs({"geom": 4326}) + assert with_crs.xindexes["geom"].crs.equals(4326) + + +def test_set_crs_arg_err(multi_geom_dataset): + with pytest.raises(ValueError, match="Cannot specify both keyword"): + multi_geom_dataset.xvec.set_crs(geom=4326, variable_crs={"geom_z": 4326}) + + +def test_set_crs_mismatch(geom_dataset): + with pytest.raises(ValueError, match="The index 'geom' already has a CRS"): + geom_dataset.xvec.set_crs(geom=4326) + + +def test_set_crs_override(geom_dataset, geom_array): + with_crs = geom_dataset.xvec.set_crs(allow_override=True, geom=4326) + # new CRS + assert with_crs.xindexes["geom"].crs.equals(4326) + # same geometries + assert shapely.equals_exact( + geom_array, + with_crs.geom, + ).all() + + +def test_set_crs_multiple(multi_geom_dataset): + with_crs = multi_geom_dataset.xvec.to_crs(geom=4326, geom_z=4326) + assert with_crs.xindexes["geom"].crs.equals(4326) + assert with_crs.xindexes["geom_z"].crs.equals(4326) diff --git a/xvec/tests/test_index.py b/xvec/tests/test_index.py index efd9fdd..53fb8f4 100644 --- a/xvec/tests/test_index.py +++ b/xvec/tests/test_index.py @@ -7,40 +7,6 @@ from xvec import GeometryIndex -@pytest.fixture(scope="session") -def geom_array(): - return np.array([shapely.Point(1, 2), shapely.Point(3, 4)]) - - -@pytest.fixture(scope="session") -def geom_dataset_no_index(geom_array): - # a dataset with a geometry coordinate but no index - ds = xr.Dataset(coords={"geom": geom_array}) - return ds.drop_indexes("geom") - - -@pytest.fixture(scope="session") -def geom_dataset(geom_dataset_no_index): - # a dataset with a geometry coordinate baked by a GeometryIndex - crs = CRS.from_user_input(26915) - return geom_dataset_no_index.set_xindex("geom", GeometryIndex, crs=crs) - - -@pytest.fixture(scope="session") -def geom_dataset_no_crs(geom_dataset_no_index): - # a dataset with a geometry coordinate baked by a GeometryIndex (no CRS) - return geom_dataset_no_index.set_xindex("geom", GeometryIndex) - - -@pytest.fixture(scope="session") -def first_geom_dataset(geom_dataset, geom_array): - return ( - xr.Dataset(coords={"geom": [geom_array[0]]}) - .drop_indexes("geom") - .set_xindex("geom", GeometryIndex, crs=geom_dataset.xindexes["geom"].crs) - ) - - def test_set_index(geom_dataset_no_index): crs = CRS.from_user_input(26915) ds = geom_dataset_no_index.set_xindex("geom", GeometryIndex, crs=crs) From 67a653576f724475b003a95b7cf486e63a7f2874 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Sat, 3 Dec 2022 17:10:26 +0100 Subject: [PATCH 07/11] remove match due to change in py3.11 --- xvec/tests/test_accessor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xvec/tests/test_accessor.py b/xvec/tests/test_accessor.py index 76d95a2..ac7726b 100644 --- a/xvec/tests/test_accessor.py +++ b/xvec/tests/test_accessor.py @@ -50,7 +50,7 @@ def test_geom_coords_names(multi_geom_dataset): assert multi_geom_dataset.xvec._geom_coords_names == ["geom", "geom_z"] assert multi_geom_dataset.xvec.geom_coords_names == ["geom", "geom_z"] - with pytest.raises(AttributeError, match="can't set attribute 'geom_coords_names'"): + with pytest.raises(AttributeError): multi_geom_dataset.xvec.geom_coords_names = ["a", "b"] From 0bfaa454eed7871ecdd63700e66a23d2b56d3406 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Sat, 3 Dec 2022 17:12:59 +0100 Subject: [PATCH 08/11] add proj to CI from conda-forge --- ci/dev.yaml | 2 ++ ci/latest.yaml | 1 + 2 files changed, 3 insertions(+) diff --git a/ci/dev.yaml b/ci/dev.yaml index 56ffe61..fd01e26 100644 --- a/ci/dev.yaml +++ b/ci/dev.yaml @@ -7,6 +7,7 @@ dependencies: # to build from source - cython - geos + - proj # testing - pytest - pytest-cov @@ -16,3 +17,4 @@ dependencies: - pip: - git+https://github.com/shapely/shapely.git@main - git+https://github.com/pydata/xarray.git@main + - git+https://github.com/pyproj4/pyproj.git diff --git a/ci/latest.yaml b/ci/latest.yaml index 9f113f4..e57ad60 100644 --- a/ci/latest.yaml +++ b/ci/latest.yaml @@ -7,6 +7,7 @@ dependencies: # required - shapely=2 - xarray + - pyproj # testing - pytest - pytest-cov From 823425f60367c3c3f83f6a47d61fe87848b2161c Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Sat, 3 Dec 2022 19:38:52 +0100 Subject: [PATCH 09/11] docs --- doc/source/api.rst | 61 ++++++++++++ doc/source/conf.py | 16 ++- environment.yml | 3 + xvec/accessor.py | 241 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 320 insertions(+), 1 deletion(-) diff --git a/doc/source/api.rst b/doc/source/api.rst index 0d27ae4..66bff71 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -4,6 +4,9 @@ API reference ============= +The API reference provides an overview of all public objects, functions and +methods and Xarray accessors implemented in Xvec. + Indexing -------- @@ -15,3 +18,61 @@ Indexing GeometryIndex GeometryIndex.crs GeometryIndex.sindex + + +.. currentmodule:: xarray + +Dataset.xvec +-------------- + +.. _dsattr: + +Properties +~~~~~~~~~~ + +.. autosummary:: + :toctree: generated/ + :template: autosummary/accessor_attribute.rst + + Dataset.xvec.geom_coords_names + + +.. _dsmeth: + +Methods +~~~~~~~ + +.. autosummary:: + :toctree: generated/ + :template: autosummary/accessor_method.rst + + Dataset.xvec.set_crs + Dataset.xvec.to_crs + + +DataArray.xvec +-------------- + +.. _daattr: + +Properties +~~~~~~~~~~ + +.. autosummary:: + :toctree: generated/ + :template: autosummary/accessor_attribute.rst + + DataArray.xvec.geom_coords_names + + +.. _dameth: + +Methods +~~~~~~~ + +.. autosummary:: + :toctree: generated/ + :template: autosummary/accessor_method.rst + + DataArray.xvec.to_crs + DataArray.xvec.set_crs \ No newline at end of file diff --git a/doc/source/conf.py b/doc/source/conf.py index a49b0da..380a47a 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -5,6 +5,15 @@ # -- Project information ----------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information +import os +import sys + +import sphinx_autosummary_accessors + + +sys.path.insert(0, os.path.abspath("../xvec/")) + +import xvec # noqa project = "Xvec" copyright = "2022, Xvec developers" @@ -20,9 +29,13 @@ "sphinx.ext.intersphinx", "myst_nb", "sphinx_copybutton", + "sphinx_autosummary_accessors", ] -templates_path = ["_templates"] +templates_path = [ + "_templates", + sphinx_autosummary_accessors.templates_path, +] exclude_patterns = [] intersphinx_mapping = { @@ -42,3 +55,4 @@ "github_url": "https://github.com/martinfleis/xvec", } nb_execution_mode = "off" +autodoc_typehints = "none" diff --git a/environment.yml b/environment.yml index 6968c44..bd8ddf3 100644 --- a/environment.yml +++ b/environment.yml @@ -22,3 +22,6 @@ dependencies: - flake8 - black - isort + - pip + - pip: + - sphinx_autosummary_accessors diff --git a/xvec/accessor.py b/xvec/accessor.py index 3e09c08..d24199d 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -11,7 +11,13 @@ @xr.register_dataarray_accessor("xvec") @xr.register_dataset_accessor("xvec") class XvecAccessor: + """Access geometry-based methods for DataArrays and Datasets with Shapely geometry. + + Currently works on coordinates with :class:`xvec.GeometryIndex`. + """ + def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): + """xvec init, nothing to be done here.""" self._obj = xarray_obj self._geom_coords_names = [ name @@ -21,6 +27,42 @@ def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): @property def geom_coords_names(self) -> List: + """Returns a list of coordinates using :class:`~xvec.GeometryIndex`. + + Returns + ------- + List + list of strings representing names of coordinates + + Examples + -------- + >>> ds = ( + ... xr.Dataset( + ... coords={ + ... "geom": np.array([shapely.Point(1, 2), shapely.Point(3, 4)]), + ... "geom_z": np.array( + ... [shapely.Point(10, 20, 30), shapely.Point(30, 40, 50)] + ... ), + ... } + ... ) + ... .drop_indexes(["geom", "geom_z"]) + ... .set_xindex("geom", xvec.GeometryIndex, crs=26915) + ... .set_xindex("geom_z", xvec.GeometryIndex, crs=26915) + ... ) + >>> ds + + Dimensions: (geom: 2, geom_z: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + * geom_z (geom_z) object POINT Z (10 20 30) POINT Z (30 40 50) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:26915) + geom_z GeometryIndex (crs=EPSG:26915) + >>> ds.xvec.geom_coords_names + ['geom', 'geom_z'] + """ return self._geom_coords_names def to_crs( @@ -28,6 +70,112 @@ def to_crs( variable_crs: Mapping[Any, Any] | None = None, **variable_crs_kwargs: Any, ): + """ + Transform :class:`shapely.Geometry` objects of a variable to a new coordinate + reference system. + + Returns a new object with all the original data in addition to the transformed + variable. The CRS the current array must be set using + :class:`~xvec.GeometryIndex`. + + This method will transform all points in all objects. It has no notion or + projecting entire geometries. All segments joining points are assumed to be + lines in the current projection, not geodesics. Objects crossing the dateline + (or other projection boundary) will have undesirable behavior. + + Parameters + ---------- + variable_crs : dict-like or None, optional + A dict where the keys are the names of the coordinates and values target + CRS in any format accepted by + :meth:`pyproj.CRS.from_user_input() ` such + as an authority string (e.g. ``"EPSG:4326"``), EPSG code (e.g. ``4326``) or + a WKT string. + **variable_crs_kwargs : optional + The keyword arguments form of ``variable_crs``. + One of ``variable_crs`` or ``variable_crs_kwargs`` must be provided. + + Returns + ------- + assigned : same type as caller + A new object with the variables transformed to target CRSs. + + See also + -------- + set_crs + + Examples + -------- + Transform coordinates backed by :class:`~xvec.GeometryIndex` from `EPSG:4326` + to `ESPG:3857`. + + >>> da = ( + ... xr.DataArray( + ... np.random.rand(2), + ... coords={"geom": [shapely.Point(1, 2), shapely.Point(3, 4)]}, + ... dims="geom", + ... ) + ... .drop_indexes("geom") + ... .set_xindex("geom", xvec.GeometryIndex, crs=4326) + ... ) + >>> da + + array([0.47575118, 0.09271935]) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + Indexes: + geom GeometryIndex (crs=EPSG:4326) + >>> da.xvec.to_crs(geom=3857) + + array([0.47575118, 0.09271935]) + Coordinates: + * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... + Indexes: + geom GeometryIndex (crs=EPSG:3857) + + The same can be done using dictionary arguments. + + >>> da.xvec.to_crs({"geom": 3857}) + + array([0.47575118, 0.09271935]) + Coordinates: + * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... + Indexes: + geom GeometryIndex (crs=EPSG:3857) + + The same applies to a :class:`xarray.Dataset`. + + >>> ds = ( + ... xr.Dataset(coords={"geom": [shapely.Point(1, 2), shapely.Point(3, 4)]}) + ... .drop_indexes("geom") + ... .set_xindex("geom", xvec.GeometryIndex, crs=4326) + ... ) + >>> ds + + Dimensions: (geom: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:4326) + >>> ds.xvec.to_crs(geom=3857) + + Dimensions: (geom: 2) + Coordinates: + * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:3857) + + Notes + ----- + Currently supports only :class:`xarray.Variable` objects that are set as + coordinates with :class:`~xvec.GeometryIndex` assigned. The implementation + currently wraps :meth:`Dataset.assign_coords ` + or :meth:`DataArray.assign_coords `. + """ if variable_crs and variable_crs_kwargs: raise ValueError( "Cannot specify both keyword and positional arguments to " @@ -96,6 +244,99 @@ def set_crs( allow_override=False, **variable_crs_kwargs: Any, ): + """Set the Coordinate Reference System (CRS) of coordinates backed by + :class:`~xvec.GeometryIndex`. + + Parameters + ---------- + variable_crs : dict-like or None, optional + A dict where the keys are the names of the coordinates and values target + CRS in any format accepted by + :meth:`pyproj.CRS.from_user_input() ` such + as an authority string (e.g. ``"EPSG:4326"``), EPSG code (e.g. ``4326``) or + a WKT string. + allow_override : bool, default False + If the the :class:`~xvec.GeometryIndex` already has a CRS, + allow to replace the existing CRS, even when both are not equal. + **variable_crs_kwargs : optional + The keyword arguments form of ``variable_crs``. + One of ``variable_crs`` or ``variable_crs_kwargs`` must be provided. + + Returns + ------- + assigned : same type as caller + A new object with the assigned target CRS. + + See also + -------- + to_crs + + Examples + -------- + The method is the most useful, when there is no CRS assigned (illustrated on + a :class:`xarray.Dataset` but the same is applicable on + a :class:`xarray.DataArray`). + + >>> ds = ( + ... xr.Dataset(coords={"geom": [shapely.Point(1, 2), shapely.Point(3, 4)]}) + ... .drop_indexes("geom") + ... .set_xindex("geom", xvec.GeometryIndex) + ... ) + >>> ds + + Dimensions: (geom: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=None) + >>> ds.xvec.set_crs(geom=4326) + + Dimensions: (geom: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:4326) + + It can also be used to overwrite the existing CRS. Note, that in most cases + you probably want to use the :meth:`to_crs` instead is such a case. + + >>> ds = ( + ... xr.Dataset(coords={"geom": [shapely.Point(1, 2), shapely.Point(3, 4)]}) + ... .drop_indexes("geom") + ... .set_xindex("geom", xvec.GeometryIndex, crs=4326) + ... ) + >>> ds + + Dimensions: (geom: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:4326) + >>> ds.xvec.set_crs(geom=3857, allow_override=True) + + Dimensions: (geom: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:3857) + + See that only the CRS has changed, not the geometries. + + + Notes + ----- + The underlying geometries are not transformed to this CRS. To + transform the geometries to a new CRS, use the :meth:`to_crs` + method. + """ if variable_crs and variable_crs_kwargs: raise ValueError( From c25839fd53516d9c9d6794ef98fa698a21d8a493 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Mon, 5 Dec 2022 14:23:18 +0100 Subject: [PATCH 10/11] Update doc/source/api.rst Co-authored-by: Benoit Bovy --- doc/source/api.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/api.rst b/doc/source/api.rst index 66bff71..8e26662 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -23,7 +23,7 @@ Indexing .. currentmodule:: xarray Dataset.xvec --------------- +------------ .. _dsattr: From b8ef560fc2833c97a026029912dae4dd2d3f5f37 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Mon, 5 Dec 2022 15:51:19 +0100 Subject: [PATCH 11/11] use is_geom_variable, geom_coords and geom_coords_indexed --- doc/source/api.rst | 8 +- xvec/accessor.py | 171 +++++++++++++++++++++++++++++++----- xvec/tests/conftest.py | 31 +++++++ xvec/tests/test_accessor.py | 53 +++++++++-- 4 files changed, 234 insertions(+), 29 deletions(-) diff --git a/doc/source/api.rst b/doc/source/api.rst index 8e26662..a16af19 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -34,7 +34,8 @@ Properties :toctree: generated/ :template: autosummary/accessor_attribute.rst - Dataset.xvec.geom_coords_names + Dataset.xvec.geom_coords + Dataset.xvec.geom_coords_indexed .. _dsmeth: @@ -46,6 +47,7 @@ Methods :toctree: generated/ :template: autosummary/accessor_method.rst + Dataset.xvec.is_geom_variable Dataset.xvec.set_crs Dataset.xvec.to_crs @@ -62,7 +64,8 @@ Properties :toctree: generated/ :template: autosummary/accessor_attribute.rst - DataArray.xvec.geom_coords_names + DataArray.xvec.geom_coords + DataArray.xvec.geom_coords_indexed .. _dameth: @@ -74,5 +77,6 @@ Methods :toctree: generated/ :template: autosummary/accessor_method.rst + DataArray.xvec.is_geom_variable DataArray.xvec.to_crs DataArray.xvec.set_crs \ No newline at end of file diff --git a/xvec/accessor.py b/xvec/accessor.py index d24199d..43f18d5 100644 --- a/xvec/accessor.py +++ b/xvec/accessor.py @@ -1,4 +1,4 @@ -from typing import Any, List, Mapping, Union +from typing import Any, Hashable, Mapping, Union import numpy as np import shapely @@ -19,20 +19,86 @@ class XvecAccessor: def __init__(self, xarray_obj: Union[xr.Dataset, xr.DataArray]): """xvec init, nothing to be done here.""" self._obj = xarray_obj - self._geom_coords_names = [ + self._geom_coords_all = [ name - for name, index in self._obj.xindexes.items() - if isinstance(index, GeometryIndex) + for name in self._obj.coords + if self.is_geom_variable(name, has_index=False) + ] + self._geom_indexes = [ + name + for name in self._obj.coords + if self.is_geom_variable(name, has_index=True) ] - @property - def geom_coords_names(self) -> List: - """Returns a list of coordinates using :class:`~xvec.GeometryIndex`. + def is_geom_variable(self, name: Hashable, has_index: bool = True): + """Check if coordinate variable is composed of :class:`shapely.Geometry`. + + Can return all such variables or only those using :class:`~xvec.GeometryIndex`. + + Parameters + ---------- + name : coordinate variable name + has_index : bool, optional + control if all variables are returned (``False``) or only those that have + :class:`xvec.GeometryIndex` assigned (``True``). By default ``False`` Returns ------- - List - list of strings representing names of coordinates + bool + + Examples + -------- + + >>> ds = ( + ... xr.Dataset( + ... coords={ + ... "geom": np.array([shapely.Point(1, 2), shapely.Point(3, 4)]), + ... "geom2": np.array([shapely.Point(1, 2), shapely.Point(3, 4)]), + ... "foo": np.array([1, 2]), + ... } + ... ) + ... .drop_indexes(["geom"]) + ... .set_xindex("geom", xvec.GeometryIndex, crs=26915) + ... ) + >>> ds + + Dimensions: (geom: 2, geom2: 2, foo: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + * geom2 (geom2) object POINT (1 2) POINT (3 4) + * foo (foo) int64 1 2 + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:26915) + >>> ds.xvec.is_geom_variable("geom") + True + >>> ds.xvec.is_geom_variable("foo") + False + >>> ds.xvec.is_geom_variable("geom2") + False + >>> ds.xvec.is_geom_variable("geom2", has_index=False) + True + + See also + -------- + geom_coords + geom_coords_indexed + """ + if isinstance(self._obj.xindexes.get(name), GeometryIndex): + return True + if not has_index: + if self._obj[name].dtype is np.dtype("O"): + # try first on a small subset + subset = self._obj[name].data[0:10] + if np.all(shapely.is_valid_input(subset)): + return np.all(shapely.is_valid_input(self._obj[name].data)) + return False + + @property + def geom_coords(self) -> Mapping[Hashable, xr.DataArray]: + """Returns a dictionary of xarray.DataArray objects corresponding to + coordinate variables composed of :class:`shapely.Geometry` objects. Examples -------- @@ -43,6 +109,9 @@ def geom_coords_names(self) -> List: ... "geom_z": np.array( ... [shapely.Point(10, 20, 30), shapely.Point(30, 40, 50)] ... ), + ... "geom_no_index": np.array( + ... [shapely.Point(1, 2), shapely.Point(3, 4)] + ... ), ... } ... ) ... .drop_indexes(["geom", "geom_z"]) @@ -51,19 +120,77 @@ def geom_coords_names(self) -> List: ... ) >>> ds - Dimensions: (geom: 2, geom_z: 2) + Dimensions: (geom: 2, geom_z: 2, geom_no_index: 2) Coordinates: - * geom (geom) object POINT (1 2) POINT (3 4) - * geom_z (geom_z) object POINT Z (10 20 30) POINT Z (30 40 50) + * geom (geom) object POINT (1 2) POINT (3 4) + * geom_z (geom_z) object POINT Z (10 20 30) POINT Z (30 40 50) + * geom_no_index (geom_no_index) object POINT (1 2) POINT (3 4) Data variables: *empty* Indexes: - geom GeometryIndex (crs=EPSG:26915) - geom_z GeometryIndex (crs=EPSG:26915) - >>> ds.xvec.geom_coords_names - ['geom', 'geom_z'] + geom GeometryIndex (crs=EPSG:26915) + geom_z GeometryIndex (crs=EPSG:26915) + >>> ds.xvec.geom_coords + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + * geom_z (geom_z) object POINT Z (10 20 30) POINT Z (30 40 50) + * geom_no_index (geom_no_index) object POINT (1 2) POINT (3 4) + + See also + -------- + geom_coords_indexed + is_geom_variable """ - return self._geom_coords_names + return self._obj[self._geom_coords_all].coords + + @property + def geom_coords_indexed(self) -> Mapping[Hashable, xr.DataArray]: + """Returns a dictionary of xarray.DataArray objects corresponding to + coordinate variables using :class:`~xvec.GeometryIndex`. + + Examples + -------- + >>> ds = ( + ... xr.Dataset( + ... coords={ + ... "geom": np.array([shapely.Point(1, 2), shapely.Point(3, 4)]), + ... "geom_z": np.array( + ... [shapely.Point(10, 20, 30), shapely.Point(30, 40, 50)] + ... ), + ... "geom_no_index": np.array( + ... [shapely.Point(1, 2), shapely.Point(3, 4)] + ... ), + ... } + ... ) + ... .drop_indexes(["geom", "geom_z"]) + ... .set_xindex("geom", xvec.GeometryIndex, crs=26915) + ... .set_xindex("geom_z", xvec.GeometryIndex, crs=26915) + ... ) + >>> ds + + Dimensions: (geom: 2, geom_z: 2, geom_no_index: 2) + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + * geom_z (geom_z) object POINT Z (10 20 30) POINT Z (30 40 50) + * geom_no_index (geom_no_index) object POINT (1 2) POINT (3 4) + Data variables: + *empty* + Indexes: + geom GeometryIndex (crs=EPSG:26915) + geom_z GeometryIndex (crs=EPSG:26915) + >>> ds.xvec.geom_coords_indexed + >>> ds.xvec.geom_coords_indexed + Coordinates: + * geom (geom) object POINT (1 2) POINT (3 4) + * geom_z (geom_z) object POINT Z (10 20 30) POINT Z (30 40 50) + + See also + -------- + geom_coords + is_geom_variable + + """ + return self._obj[self._geom_indexes].coords def to_crs( self, @@ -122,14 +249,14 @@ def to_crs( array([0.47575118, 0.09271935]) Coordinates: - * geom (geom) object POINT (1 2) POINT (3 4) + * geom (geom) object POINT (1 2) POINT (3 4) Indexes: geom GeometryIndex (crs=EPSG:4326) >>> da.xvec.to_crs(geom=3857) array([0.47575118, 0.09271935]) Coordinates: - * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... + * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... Indexes: geom GeometryIndex (crs=EPSG:3857) @@ -163,7 +290,7 @@ def to_crs( Dimensions: (geom: 2) Coordinates: - * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... + * geom (geom) object POINT (111319.49079327357 222684.20850554405) POIN... Data variables: *empty* Indexes: @@ -286,7 +413,7 @@ def set_crs( Dimensions: (geom: 2) Coordinates: - * geom (geom) object POINT (1 2) POINT (3 4) + * geom (geom) object POINT (1 2) POINT (3 4) Data variables: *empty* Indexes: @@ -295,7 +422,7 @@ def set_crs( Dimensions: (geom: 2) Coordinates: - * geom (geom) object POINT (1 2) POINT (3 4) + * geom (geom) object POINT (1 2) POINT (3 4) Data variables: *empty* Indexes: diff --git a/xvec/tests/conftest.py b/xvec/tests/conftest.py index c67b9e8..c3e6455 100644 --- a/xvec/tests/conftest.py +++ b/xvec/tests/conftest.py @@ -59,3 +59,34 @@ def multi_geom_dataset(geom_array, geom_array_z): .set_xindex("geom", GeometryIndex, crs=26915) .set_xindex("geom_z", GeometryIndex, crs=26915) ) + + +@pytest.fixture(scope="session") +def multi_geom_no_index_dataset(geom_array, geom_array_z): + return ( + xr.Dataset( + coords={ + "geom": geom_array, + "geom_z": geom_array_z, + "geom_no_ix": geom_array, + } + ) + .drop_indexes(["geom", "geom_z"]) + .set_xindex("geom", GeometryIndex, crs=26915) + .set_xindex("geom_z", GeometryIndex, crs=26915) + ) + + +@pytest.fixture(scope="session") +def multi_geom_one_ix_foo(geom_array): + return ( + xr.Dataset( + coords={ + "geom": geom_array, + "geom2": geom_array, + "foo": np.array([1, 2]), + } + ) + .drop_indexes(["geom"]) + .set_xindex("geom", GeometryIndex, crs=26915) + ) diff --git a/xvec/tests/test_accessor.py b/xvec/tests/test_accessor.py index ac7726b..14e3b86 100644 --- a/xvec/tests/test_accessor.py +++ b/xvec/tests/test_accessor.py @@ -43,15 +43,58 @@ def test_accessor(multi_geom_dataset): xr.testing.assert_identical(multi_geom_dataset.xvec._obj, multi_geom_dataset) -# Test .xvec geom_coords_names +# Test .xvec geom_coords -def test_geom_coords_names(multi_geom_dataset): - assert multi_geom_dataset.xvec._geom_coords_names == ["geom", "geom_z"] - assert multi_geom_dataset.xvec.geom_coords_names == ["geom", "geom_z"] +def test_geom_coords(multi_geom_no_index_dataset): + assert multi_geom_no_index_dataset.xvec._geom_coords_all == [ + "geom", + "geom_z", + "geom_no_ix", + ] + actual = multi_geom_no_index_dataset.xvec.geom_coords + expected = multi_geom_no_index_dataset.coords + actual.keys() == expected.keys() + + # check assignment with pytest.raises(AttributeError): - multi_geom_dataset.xvec.geom_coords_names = ["a", "b"] + multi_geom_no_index_dataset.xvec.geom_coords = ( + multi_geom_no_index_dataset.coords + ) + + +def test_geom_coords_indexed(multi_geom_no_index_dataset): + assert multi_geom_no_index_dataset.xvec._geom_indexes == ["geom", "geom_z"] + + actual = multi_geom_no_index_dataset.xvec.geom_coords_indexed + expected = multi_geom_no_index_dataset.coords + actual.keys() == expected.keys() + + # check assignment + with pytest.raises(AttributeError): + multi_geom_no_index_dataset.xvec.geom_coords = ( + multi_geom_no_index_dataset.coords + ) + + +# Test .xvec.is_geom_variable +@pytest.mark.parametrize( + "label,has_index,expected", + [ + ("geom", True, True), + ("geom", False, True), + ("geom2", True, False), + ("geom2", False, True), + ("foo", True, False), + ("foo", False, False), + ], +) +def test_is_geom_variable(multi_geom_one_ix_foo, label, has_index, expected): + assert ( + multi_geom_one_ix_foo.xvec.is_geom_variable(label, has_index=has_index) + == expected + ) # Test .xvec.to_crs