Repository: earth-chris/elapid

Issue:
Report progress for raster annotation
The `raster_values_from_geoseries()` function has iterables for reading data 1) for each raster input and 2) for each point feature. We should add some way to track the progress for each of these with `tqdm`.

Since one of the iterations is via `df.apply()`, we could use [this feature](https://stackoverflow.com/questions/18603270/progress-indicator-during-pandas-operations).

Discussion:

Patch:
diff --git a/elapid/geo.py b/elapid/geo.py
index 7a9481e..6b3111f 100644
--- a/elapid/geo.py
+++ b/elapid/geo.py
@@ -20,6 +20,7 @@ from elapid.utils import (
 )
 
 tqdm = get_tqdm()
+tqdm.pandas(desc="Geometry", leave=False)
 
 
 def xy_to_geoseries(x, y, crs="epsg:4236"):
@@ -33,6 +34,11 @@ def xy_to_geoseries(x, y, crs="epsg:4236"):
     Returns:
         gs: Point geometry geoseries
     """
+    if not hasattr(x, "__iter__"):
+        x = [x]
+    if not hasattr(y, "__iter__"):
+        y = [y]
+
     points = [Point(x, y) for x, y in zip(x, y)]
     gs = gpd.GeoSeries(points, crs=crs)
 
@@ -151,6 +157,71 @@ def pseudoabsence_from_geoseries(geoseries, count, overestimate=2):
     return points
 
 
+def parse_crs_string(string):
+    """Parses a string to determine the CRS/spatial projection format.
+
+    Args:
+        string: a string with CRS/projection data.
+
+    Returns:
+        crs_type: Str in ["wkt", "proj4", "epsg", "string"].
+    """
+    if "epsg:" in string.lower():
+        return "epsg"
+    elif "+proj" in string:
+        return "proj4"
+    elif "SPHEROID" in string:
+        return "wkt"
+    else:
+        return "string"
+
+
+def string_to_crs(string):
+    """Converts a crs/projection string to a pyproj-readable CRS object
+
+    Args:
+        string: a crs/projection string.
+        crs_type: the type of crs/projection string, in ["wkt", "proj4", "epsg", "string"].
+
+    Returns:
+        crs: a rasterio.crs.CRS object
+    """
+    crs_type = parse_crs_string(string)
+
+    if crs_type == "epsg":
+        auth, code = string.split(":")
+        crs = rio.crs.CRS.from_epsg(int(code))
+    elif crs_type == "proj4":
+        crs = rio.crs.CRS.from_proj4(string)
+    elif crs_type == "wkt":
+        crs = rio.crs.CRS.from_wkt(string)
+    else:
+        crs = rio.crs.CRS.from_string(string)
+
+    return crs
+
+
+def crs_match(crs1, crs2):
+    """Determines whether two coordinate reference systems are the same.
+
+    Args:
+        crs1: the first CRS, from a rasterio dataset, a GeoDataFrame, or a string with projection parameters.
+        crs2: the second CRS, from the same sources above.
+
+    Returns:
+        matches: Boolean for whether the CRS match.
+    """
+    # normalize string inputs via rasterio
+    if type(crs1) is str:
+        crs1 = string_to_crs(crs1)
+    if type(crs2) is str:
+        crs2 = string_to_crs(crs2)
+
+    matches = crs1 == crs2
+
+    return matches
+
+
 def raster_values_from_vector(vector_path, raster_paths, labels=None, drop_na=True):
     """Reads and stores pixel values from rasters using a point-format vector file.
 
@@ -203,7 +274,7 @@ def raster_values_from_geoseries(geoseries, raster_paths, labels=None, drop_na=T
 
     # apply this function over every geoseries row
     def read_pixel_value(point, source):
-        row, col = source.index(point.x, point.y)
+        row, col = source.index(point.geometry.x, point.geometry.y)
         window = rio.windows.Window(col, row, 1, 1)
         value = source.read(window=window)
         return np.squeeze(value)
@@ -211,8 +282,14 @@ def raster_values_from_geoseries(geoseries, raster_paths, labels=None, drop_na=T
     # apply this function over every raster_path
     def parallel_raster_reads(raster_path):
         with rio.open(raster_path) as src:
-            points = geoseries.to_crs(src.crs)
-            values = points.apply(read_pixel_value, source=src)
+
+            # reproject if necessary
+            if not crs_match(geoseries.crs, src.crs):
+                points = geoseries.to_crs(src.crs).to_frame("geometry")
+            else:
+                points = geoseries.to_frame("geometry")
+
+            values = points.progress_apply(read_pixel_value, axis=1, result_type="expand", source=src)
             if drop_na and src.nodata is not None:
                 values.replace(src.nodata, np.NaN, inplace=True)
 
diff --git a/environment.yml b/environment.yml
index 9711b6f..b6c154d 100644
--- a/environment.yml
+++ b/environment.yml
@@ -5,7 +5,7 @@ channels:
   - defaults
 
 dependencies:
-  - python==3.7
+  - python>=3.7
   - descartes
   - gdal
   - geopandas
@@ -15,7 +15,6 @@ dependencies:
   - rasterio>=1.0
   - rtree
   - scikit-learn
-  - scipy
   - shapely
   - tqdm
   - pip
@@ -24,3 +23,4 @@ dependencies:
     - mkdocs-material
     - mkdocs-jupyter
     - mkdocstrings
+    - livereload


Test patch:
diff --git a/tests/test_geo.py b/tests/test_geo.py
new file mode 100644
index 0000000..57ac775
--- /dev/null
+++ b/tests/test_geo.py
@@ -0,0 +1,120 @@
+import os
+from copy import copy
+
+import geopandas as gpd
+import numpy as np
+import rasterio as rio
+
+from elapid import geo
+
+# set the test raster data paths
+directory_path, script_path = os.path.split(os.path.abspath(__file__))
+data_path = os.path.join(directory_path, "data")
+raster_1b = os.path.join(data_path, "test-raster-1band.tif")
+raster_2b = os.path.join(data_path, "test-raster-2bands.tif")
+raster_1b_offset = os.path.join(data_path, "test-raster-1band-offset.tif")
+with rio.open(raster_1b, "r") as src:
+    raster_1b_profile = copy(src.profile)
+
+# set some crs variables
+wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
+proj4 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
+epsg_code = 4326
+epsg_str = f"EPSG:{epsg_code}"
+
+
+def test_xy_to_geoseries():
+    lon, lat = (-120, 38)
+    geoseries = geo.xy_to_geoseries(lon, lat, crs=epsg_str)
+    assert geoseries.x[0] == lon
+    assert geoseries.y[0] == lat
+
+
+# TODO
+def test_pseudoabsence_from_raster():
+    pass
+
+
+# TODO
+def test_pseudoabsence_from_bias_file():
+    pass
+
+
+# TODO
+def test_pseudoabsence_from_vector():
+    pass
+
+
+# TODO
+def test_pseudoabsence_from_geoseries():
+    pass
+
+
+def test_parse_crs_string():
+    assert geo.parse_crs_string(wkt) == "wkt"
+    assert geo.parse_crs_string(proj4) == "proj4"
+    assert geo.parse_crs_string(epsg_str) == "epsg"
+
+
+def test_string_to_crs():
+    crs = geo.string_to_crs(wkt)
+    assert crs.to_epsg() == epsg_code
+
+    crs = geo.string_to_crs(proj4)
+    assert crs.to_epsg() == epsg_code
+
+    crs = geo.string_to_crs(epsg_str)
+    assert crs.to_epsg() == epsg_code
+
+
+def test_crs_match():
+    lon, lat = (-120, 38)
+    geoseries = geo.xy_to_geoseries(lon, lat, crs=epsg_str)
+    gpd_crs = geoseries.crs
+    rio_crs = rio.crs.CRS.from_epsg(epsg_code)
+
+    assert geo.crs_match(rio_crs, gpd_crs) is True
+    assert geo.crs_match(gpd_crs, wkt) is True
+    assert geo.crs_match(rio_crs, epsg_str) is True
+    assert geo.crs_match(wkt, proj4) is True
+
+    assert geo.crs_match(rio_crs, "epsg:32610") is False
+
+
+# TODO
+def test_raster_values_from_vector():
+    pass
+
+
+def test_raster_values_from_geoseries():
+    # create a single point in the origin of the test data
+    with rio.open(raster_1b, "r") as src:
+        x, y = src.xy(0, 0)
+        geoseries = geo.xy_to_geoseries(x, y, crs=src.crs)
+
+    # test on one band input
+    df = geo.raster_values_from_geoseries(geoseries, [raster_1b])
+    b1 = df["band_1"].iloc[0]
+    assert b1 == 0
+
+    # test on two band input
+    df = geo.raster_values_from_geoseries(geoseries, [raster_2b])
+    b1 = df["band_1"].iloc[0]
+    b2 = df["band_2"].iloc[0]
+    assert b1 == 0
+    assert b2 == 65280
+
+    # test on multi-raster input
+    df = geo.raster_values_from_geoseries(geoseries, [raster_1b, raster_2b])
+    b1 = df["band_1"].iloc[0]
+    b2 = df["band_2"].iloc[0]
+    b3 = df["band_3"].iloc[0]
+    assert b1 == 0
+    assert b2 == 0
+    assert b3 == 65280
+
+    # test it works on lots of points
+    n = 100
+    df = geo.raster_values_from_geoseries(geoseries.repeat(n), [raster_1b])
+    b1 = df["band_1"].iloc[n - 1]
+    assert b1 == 0