ok

ecoscape_layers/constants.py

@@ -1,5 +1,5 @@
 RESAMPLING_METHODS = ["near", "bilinear", "cubic", "cubicspline", "lanczos", "average", "rms", "mode", "max", "min", "med", "q1", "q3", "sum"]
-REFINE_METHODS = ["forest", "forest_add308", "allsuitable", "majoronly"]
+REFINE_METHODS = ["forest", "forest_africa", "forest_add308", "allsuitable", "majoronly"]
 
 HAB_308 = {
     "code": "3.8",
@@ -9,4 +9,4 @@
     "majorimportance": "Yes",
     "map_code": 308,
     "resistance": 0
-}
+}

ecoscape_layers/layers.py

@@ -15,6 +15,7 @@
 
 import tqdm
 import psutil
+from itertools import chain
 
 
 class LayerGenerator(object):
@@ -23,6 +24,26 @@ class LayerGenerator(object):
     This class maintains a common CRS, resolution, and resampling method for this purpose.
     """
 
+    # IUCN map code ranges
+    # NOTE: These codes can be found on https://www.iucnredlist.org/resources/habitat-classification-scheme
+    FORESTS = (100, 200)
+    SAVANNAS = (200, 300)
+    SHRUBLANDS = (300, 400)
+    GRASSLANDS = (400, 500)
+    WETLANDS = (500, 600)
+    ROCKY = (600, 700)
+    CAVES = (700, 800)
+    DESERTS = (800, 900)
+    MARINE_NERITIC = (900, 1000)
+    MARINE_OCEANIC = (1000, 1100)
+    MARINE_DEEP = (1100, 1200)
+    MARINE_INTERTIDAL = (1200, 1300)
+    MARINE_COSTAL = (1300, 1400)
+    ARTIFICIAL = (1400, 1600)
+    INTRODUCED_VEGETATION = (1600, 1700)
+    OTHER = (1700, 1800)
+    UNKNOWN = (1800, 1900)
+
     def __init__(
         self,
         redlist_key: str,
@@ -53,53 +74,17 @@ def get_map_codes(self, landcover: GeoTiff) -> list[int]:
         :param landcover: The landcover map in GeoTiff format.
         """
 
-        # get the bounds of the landcover map
-        bounds = landcover.dataset.bounds
-
-        # convert bounds to id num
-        bounds_id = hash(bounds)
-
-        # get the file name of the landcover map
-        landcover_fn_basename = os.path.basename(self.landcover_fn)
-        landcover_fn_basename = os.path.splitext(landcover_fn_basename)[0]
-
-        # add the bounds to the landcover file name
-        map_codes_id = f"{landcover_fn_basename}_{bounds_id}"
-
-        # create map codes file name
-        map_codes_fn = f"map_codes_{map_codes_id}.csv"
-
-        # create map codes directory
-        map_codes_dir = os.path.join(os.path.dirname(self.landcover_fn), "map_codes")
-        os.makedirs(map_codes_dir, exist_ok=True)
-
-        # create the full path to the map codes file
-        map_codes_output_path = os.path.join(map_codes_dir, map_codes_fn)
-
-        # check if the map codes file exists
-        if os.path.exists(map_codes_output_path):
-            # read the map codes file
-            with open(map_codes_output_path, "r") as csvfile:
-                reader = csv.reader(csvfile)
-                next(reader)
-                map_codes = [int(row[0]) for row in reader]
-
-            return map_codes
-
         # get all unique map codes from the landcover map
         tile = landcover.get_all_as_tile()
 
         if tile is None:
             raise ValueError("Landcover file is empty.")
 
-        map_codes: list[int] = sorted(list(np.unique(tile.m)))
+        # get map data from tile.m as a numpy matrix
+        map_data: np.ndarray = tile.m
 
-        # write the map codes to a file
-        with open(map_codes_output_path, "w", newline="") as csvfile:
-            writer = csv.writer(csvfile)
-            writer.writerow(["map_code"])
-            for map_code in map_codes:
-                writer.writerow([map_code])
+        # get unique map codes
+        map_codes: list[int] = sorted(list(set(chain(map_data.flatten()))))
 
         return map_codes
 
@@ -191,7 +176,8 @@ def generate_resistance_table(
         """
 
         # get the map codes
-        map_codes = self.get_map_codes(landcover)
+        # using range 2000 as this is all available map codes
+        map_codes = range(2000)
 
         with open(output_path, "w", newline="") as csvfile:
             writer = csv.DictWriter(csvfile, fieldnames=habitats[0].keys())
@@ -202,7 +188,11 @@ def generate_resistance_table(
                 h = code_to_habitat.get(map_code)
                 if h is not None:
                     if refine_method == "forest" or refine_method == "forest_add308":
-                        h["resistance"] = 0 if map_code >= 100 and map_code < 200 else h["resistance"]
+                        h["resistance"] = 0 if in_hab(self.FORESTS) else h["resistance"]
+                    elif refine_method == "forest_africa":
+                        # will have the forest and shrublands set to 0 regardless of the habitat data
+                        h["resistance"] = 0 if in_hab(self.FORESTS, self.SHRUBLANDS) else h["resistance"]
+
                     writer.writerow(h.values())
                 else:
                     default_row = {'map_code': map_code}, {'resistance': 0 if map_code >= 100 and map_code < 200 else 1}
@@ -218,14 +208,28 @@ def get_good_terrain(self, habitats, refine_method="forest_add308") -> list[int]
         :return: list of map codes filtered by refine_method.
         """
 
+        # Predefine good terrain map codes
+        suit_hab = [h["map_code"] for h in habitats if h["suitability"] == "Suitable"]
+        maj_hab = [h["map_code"] for h in habitats if h["majorimportance"] == "Yes"]
+
+        # function to combine ranges from range function
+        def comb_ranges(*ranges: tuple[int, int]) -> list[int]:
+            codes = []
+            for r in ranges:
+                codes.extend(list(range(r[0], r[1])))
+            return codes
+
         if refine_method == "forest":
-            return [x for x in range(100, 110)]
+            return comb_ranges(self.FORESTS)
         elif refine_method == "forest_add308":
-            return [x for x in range(100, 110)] + [308]
+            return comb_ranges(self.FORESTS) + [308]
+        elif refine_method == "forest_africa":
+            hab_ranges = comb_ranges(self.FORESTS, self.SHRUBLANDS, self.INTRODUCED_VEGETATION)
+            return list(set(hab_ranges + suit_hab + maj_hab))
         elif refine_method == "allsuitable":
-            return [hab["map_code"] for hab in habitats if hab["suitability"] == "Suitable"]
+            return suit_hab
         elif refine_method == "majoronly":
-            return [hab["map_code"] for hab in habitats if hab["majorimportance"] == "Yes"]
+            return maj_hab
         else:
             return []
 

ecoscape_layers/redlist.py

@@ -15,7 +15,7 @@ def __init__(self, redlist_key: str, ebird_key: str | None = None):
         self.redlist_params = {"token": redlist_key}
         self.ebird_key = ebird_key
 
-    def get_from_redlist(self, url) -> list[dict[str, str | int | float]]:
+    def get_from_redlist(self, url):
         """
         Convenience function for sending GET request to Red List API with the key.
 
@@ -67,7 +67,7 @@ def get_habitat_data(
         """
         url = f"https://apiv3.iucnredlist.org/api/v3/habitats/species/name/{name}"
         if region is not None:
-            url += "/region/{1}".format(region)
+            url += f"/region/{region}"
 
         habs = self.get_from_redlist(url)
 

ecoscape_layers/utils.py

@@ -64,8 +64,8 @@ def transform_box(
     min_lat: float,
     max_lon: float,
     max_lat: float,
+    crs_in: str,
     crs_out: str,
-    crs_in: str = "EPSG:4326",
 ):
     """
     Transforms a bounding box from one coordinate reference system (CRS) to another.
@@ -75,8 +75,8 @@ def transform_box(
         min_lat (float): The minimum latitude of the bounding box.
         max_lon (float): The maximum longitude of the bounding box.
         max_lat (float): The maximum latitude of the bounding box.
-        crs_out (str): The output CRS in EPSG format.
-        crs_in (str): The input CRS in EPSG format.
+        crs_in (str): The input CRS in EPSG format. For example "EPSG:4326".
+        crs_out (str): The output CRS in EPSG format. For example "EPSG:3395".
 
     Returns:
         tuple: A tuple containing the transformed coordinates of the bounding box in the output CRS.

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "ecoscape-layers"
-version = "0.0.8"
+version = "0.0.9"
 authors = [
   {name="Jasmine Tai", email="[email protected]"},
   {name="Coen Adler", email="[email protected]"},