Add point set matching for trees

GDRT/constants.py

@@ -1,6 +1,9 @@
 import typing
 from pathlib import Path
 
+import pyproj
+
 PATH_TYPE = typing.Union[Path, str]
 
 DATA_FOLDER = Path(Path(__file__).parent, "..", "data").resolve()
+LAT_LON_CRS = pyproj.CRS.from_epsg(4326)

GDRT/geospatial_utils.py

@@ -7,6 +7,8 @@
 from contourpy import contour_generator
 from rasterio.features import rasterize
 
+from GDRT.constants import LAT_LON_CRS
+
 
 def mask_to_shapely(mask: np.ndarray) -> shapely.MultiPolygon:
     """Convert a binary mask to a shapely polygon representing the positive regions
@@ -113,3 +115,36 @@ def extract_largest_oriented_rectangle(gpd_multipolygon, raster_resolution=1.0):
     rect_gpd = gpd.GeoDataFrame(geometry=[rect_shapely], crs=gpd_multipolygon.crs)
 
     return rect_gpd
+
+
+def get_projected_CRS(lat, lon, assume_western_hem=True):
+    if assume_western_hem and lon > 0:
+        lon = -lon
+    # https://gis.stackexchange.com/questions/190198/how-to-get-appropriate-crs-for-a-position-specified-in-lat-lon-coordinates
+    epgs_code = 32700 - round((45 + lat) / 90) * 100 + round((183 + lon) / 6)
+    crs = pyproj.CRS.from_epsg(epgs_code)
+    return crs
+
+
+def ensure_projected_CRS(geodata: gpd.GeoDataFrame):
+    """Returns a projected geodataframe from the provided geodataframe by converting it to
+    ESPG:4326 (if not already) and determining the projected CRS from the point
+    coordinates.
+
+    Args:
+        geodata (gpd.GeoDataGrame): Original geodataframe that is potentially unprojected
+    Returns:
+        gpd.GeoDataGrame: projected geodataframe
+    """
+    # If CRS is projected return immediately
+    if geodata.crs.is_projected:
+        return geodata
+
+    # If CRS is geographic and not long-lat, convert it to long-lat
+    if geodata.crs.is_geographic and geodata.crs != LAT_LON_CRS:
+        geodata = geodata.to_crs(LAT_LON_CRS)
+
+    # Convert geographic long-lat CRS to projected CRS
+    point = geodata["geometry"].iloc[0].centroid
+    geometric_crs = get_projected_CRS(lon=point.x, lat=point.y)
+    return geodata.to_crs(geometric_crs)