Add Vector.reproject()

geoutils/georaster/raster.py

@@ -1446,13 +1446,11 @@ def reproject(
         """
 
         # Check that either dst_ref or dst_crs is provided
-        if dst_ref is not None:
-            if dst_crs is not None:
-                raise ValueError("Either of `dst_ref` or `dst_crs` must be set. Not both.")
-        else:
-            # In case dst_res or dst_size is set, use original CRS
-            if dst_crs is None:
-                dst_crs = self.crs
+        if dst_ref is not None and dst_crs is not None:
+            raise ValueError("Either of `dst_ref` or `dst_crs` must be set. Not both.")
+        # If none are provided, simply preserve the CRS
+        elif dst_ref is None and dst_crs is None:
+            dst_crs = self.crs
 
         # Case a raster is provided as reference
         if dst_ref is not None:
@@ -1464,11 +1462,12 @@ def reproject(
             elif isinstance(dst_ref, rio.io.MemoryFile) or isinstance(dst_ref, rasterio.io.DatasetReader):
                 ds_ref = dst_ref
             elif isinstance(dst_ref, str):
-                assert os.path.exists(dst_ref), "Reference raster does not exist"
+                if not os.path.exists(dst_ref):
+                    raise ValueError("Reference raster does not exist.")
                 ds_ref = Raster(dst_ref, load_data=False)
             else:
-                raise ValueError(
-                    "Type of dst_ref not understood, must be path to file (str), Raster or rasterio data set"
+                raise TypeError(
+                    "Type of dst_ref not understood, must be path to file (str), Raster or rasterio data set."
                 )
 
             # Read reprojecting params from ref raster

geoutils/geovector.py

@@ -3,16 +3,19 @@
 """
 from __future__ import annotations
 
+import os
 import pathlib
 import warnings
 from collections import abc
 from numbers import Number
 from typing import Literal, TypeVar, overload
 
+import fiona
 import geopandas as gpd
 import matplotlib.pyplot as plt
 import numpy as np
 import rasterio as rio
+import rasterio.errors
 import shapely
 from rasterio import features, warp
 from rasterio.crs import CRS
@@ -53,7 +56,7 @@ def __init__(self, filename: str | pathlib.Path | gpd.GeoDataFrame):
         else:
             raise ValueError("filename argument not recognised.")
 
-        self.crs = self.ds.crs
+        self._crs = self.ds.crs
 
     def __repr__(self) -> str:
         return str(self.ds.__repr__())
@@ -89,6 +92,12 @@ def bounds(self) -> rio.coords.BoundingBox:
         """Get a bounding box of the total bounds of the Vector."""
         return rio.coords.BoundingBox(*self.ds.total_bounds)
 
+    @property
+    def crs(self) -> rio.crs.CRS:
+        """Ensure CRS cannot be set"""
+        self._crs = self.ds.crs
+        return self._crs
+
     def copy(self: VectorType) -> VectorType:
         """Return a copy of the Vector."""
         # Utilise the copy method of GeoPandas
@@ -146,6 +155,62 @@ def crop(
             new_vector.ds = new_vector.ds.cx[xmin:xmax, ymin:ymax]
             return new_vector
 
+    def reproject(
+        self: Vector,
+        dst_ref: gu.Raster | rio.io.DatasetReader | VectorType | gpd.GeoDataFrame | str | None = None,
+        dst_crs: CRS | str | int | None = None,
+    ) -> Vector:
+        """
+        Reproject vector to a specified CRS and, optionally, cropped to certain bounds (no cropping by default).
+
+        The output CRS can either be given by a reference Raster or Vector (using `dst_ref`) or by manually
+        providing the output CRS (`dst_crs`)
+
+        To reproject a Vector with different source bounds, first run Vector.crop().
+
+        :param dst_ref: A reference raster or vector. If set will reproject to the CRS of this reference, and
+            optionally crop to the same bounds. Can be provided as Raster, Vector, rasterio dataset,
+            geopandas dataframe, or path to the file.
+        :param dst_crs: Specify the Coordinate Reference System or EPSG to reproject to. If dst_ref not set,
+            defaults to self.crs.
+
+        :returns: Reprojected vector
+        """
+
+        # Check that either dst_ref or dst_crs is provided
+        if (dst_ref is not None and dst_crs is not None) or (dst_ref is None and dst_crs is None):
+            raise ValueError("Either of `dst_ref` or `dst_crs` must be set. Not both.")
+
+        # Case a raster or vector is provided as reference
+        if dst_ref is not None:
+
+            # Check that dst_ref type is either str, Raster or rasterio data set
+            # Preferably use Raster instance to avoid rasterio data set to remain open. See PR #45
+            if isinstance(dst_ref, (gu.Raster, gu.Vector)):
+                ds_ref = dst_ref
+            elif isinstance(dst_ref, (rio.io.DatasetReader, gpd.GeoDataFrame)):
+                ds_ref = dst_ref
+            elif isinstance(dst_ref, str):
+                if not os.path.exists(dst_ref):
+                    raise ValueError("Reference raster or vector path does not exist.")
+                try:
+                    ds_ref = gu.Raster(dst_ref, load_data=False)
+                except rasterio.errors.RasterioIOError:
+                    try:
+                        ds_ref = Vector(dst_ref)
+                    except fiona.errors.DriverError:
+                        raise ValueError("Could not open raster or vector with rasterio or fiona.")
+            else:
+                raise TypeError("Type of dst_ref must be string path to file, Raster or Vector.")
+
+            # Read reprojecting params from ref raster
+            dst_crs = ds_ref.crs
+        else:
+            # Determine user-input target CRS
+            dst_crs = CRS.from_user_input(dst_crs)
+
+        return Vector(self.ds.to_crs(crs=dst_crs))
+
     def create_mask(
         self,
         rst: str | gu.Raster | None = None,

tests/test_geovector.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import pathlib
+import re
 
 import geopandas as gpd
 import numpy as np
@@ -75,6 +76,48 @@ def test_bounds(self) -> None:
         assert bounds.right == self.glacier_outlines.ds.total_bounds[2]
         assert bounds.top == self.glacier_outlines.ds.total_bounds[3]
 
+    def test_reproject(self) -> None:
+        """Test that the reproject function works as intended"""
+
+        v0 = gu.Vector(self.aster_outlines_path)
+        r0 = gu.Raster(self.aster_dem_path)
+        v1 = gu.Vector(self.everest_outlines_path)
+
+        # First, test with a EPSG integer
+        v1 = v0.reproject(dst_crs=32617)
+        assert isinstance(v1, gu.Vector)
+        assert v1.crs.to_epsg() == 32617
+
+        # Check that the reprojection is the same as with geopandas
+        gpd1 = v0.ds.to_crs(epsg=32617)
+        assert_geodataframe_equal(gpd1, v1.ds)
+
+        # Second, with a Raster object
+        v2 = v0.reproject(r0)
+        assert v2.crs == r0.crs
+
+        # Third, with a Vector object that has a different CRS
+        assert v0.crs != v1.crs
+        v3 = v0.reproject(v1)
+        assert v3.crs == v1.crs
+
+        # Fourth, check that errors are raised when appropriate
+        # When no destination CRS is defined, or both dst_crs and dst_ref are passed
+        with pytest.raises(ValueError, match=re.escape("Either of `dst_ref` or `dst_crs` must be set. Not both.")):
+            v0.reproject()
+            v0.reproject(dst_ref=r0, dst_crs=32617)
+        # If the path provided does not exist
+        with pytest.raises(ValueError, match=re.escape("Reference raster or vector path does not exist.")):
+            v0.reproject(dst_ref="tmp.lol")
+        # If it exists but cannot be opened by rasterio or fiona
+        with pytest.raises(ValueError, match=re.escape("Could not open raster or vector with rasterio or fiona.")):
+            v0.reproject(dst_ref="geoutils/examples.py")
+        # If input of wrong type
+        with pytest.raises(
+            TypeError, match=re.escape("Type of dst_ref must be string path to file, Raster or Vector.")
+        ):
+            v0.reproject(dst_ref=10)  # type: ignore
+
     def test_rasterize_proj(self) -> None:
 
         # Capture the warning on resolution not matching exactly bounds