Update habitat creation

config/config.yaml

@@ -0,0 +1,28 @@
+# LIFE Pipeline Configuration
+# ===========================
+
+# Taxonomic classes to process
+taxa:
+  - AMPHIBIA
+  - AVES
+  - MAMMALIA
+  - REPTILIA
+
+# Pixel scale
+pixel_scale: 0.016666666666667
+
+# Hyde projection data
+hyde_projection: '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.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
+hyde_pixel_scale: 0.08333333333333333
+
+# Zenodo configuration for downloading raw habitat
+zenodo:
+    jung_habitat:
+        zenodo_id: 4058819
+        filename: "iucn_habitatclassification_composite_lvl2_ver004.zip"
+    jung_habitat_updates:
+        zenodo_id: 4058819
+        filename: "lvl2_changemasks_ver004.zip"
+    jung_pnv:
+        zenodo_id: 4038749
+        filename: "pnv_lvl1_004.zip"

deltap/delta_p_scaled.py

@@ -4,66 +4,57 @@
 from pathlib import Path
 
 import pandas as pd
-import yirgacheffe.operators as yo
-from yirgacheffe.layers import RasterLayer
+import yirgacheffe as yg
 
 SCALE = 1e6
 
 def delta_p_scaled_area(
     input_path: Path,
     diff_area_map_path: Path,
-    totals_path: Path,
+    species_totals_path: Path,
     output_path: Path,
 ):
     os.makedirs(output_path.parent, exist_ok=True)
 
     per_taxa = [
-        RasterLayer.layer_from_file(os.path.join(input_path, x))
-        for x in sorted(input_path.glob("*.tif"))
+        yg.read_raster(x) for x in sorted(input_path.glob("*.tif"))
     ]
     if not per_taxa:
         sys.exit(f"Failed to find any per-taxa maps in {input_path}")
 
-    area_restore = RasterLayer.layer_from_file(diff_area_map_path)
+    species_total_counts = pd.read_csv(species_totals_path)
 
-    total_counts = pd.read_csv(totals_path)
+    diff_area = yg.read_raster(diff_area_map_path)
+    diff_area_rescaled = yg.where(diff_area < SCALE, float('nan'), diff_area / SCALE)
 
-    area_restore_filter = yo.where(area_restore < SCALE, float('nan'), area_restore) / SCALE
-
-    with RasterLayer.empty_raster_layer_like(
-        area_restore,
-        filename=output_path,
-        nodata=float('nan'),
-        bands=len(per_taxa) + 1
-    ) as result:
+    # Process all species in total
+    total_species_count = int(species_total_counts[species_total_counts.taxa=="all"]["count"].values[0])
+    summed_layer = yg.sum(per_taxa)
+    all_scaled_filtered_layer = yg.where(
+        diff_area_rescaled != 0,
+        ((summed_layer / diff_area_rescaled) * -1.0) / total_species_count,
+        float('nan')
+    )
 
-        species_count = int(total_counts[total_counts.taxa=="all"]["count"].values[0])
+    bands = [all_scaled_filtered_layer]
+    labels = ["all"]
 
-        result._dataset.GetRasterBand(1).SetDescription("all")  # pylint: disable=W0212
-        summed_layer = per_taxa[0]
-        for layer in per_taxa[1:]:
-            summed_layer = summed_layer + layer
+    # Now per taxa
+    for inlayer in per_taxa:
+        # get the taxa from the filename
+        _, name = os.path.split(inlayer.name)
+        taxa = name[:-4]
+        labels.append(taxa)
 
-        scaled_filtered_layer = yo.where(
-            area_restore_filter != 0,
-            ((summed_layer / area_restore_filter) * -1.0) / species_count,
+        taxa_species_count = int(species_total_counts[species_total_counts.taxa==taxa]["count"].values[0])
+        scaled_filtered_layer = yg.where(
+            diff_area_rescaled != 0,
+            ((inlayer / diff_area_rescaled) * -1.0) / taxa_species_count,
             float('nan')
         )
-        scaled_filtered_layer.parallel_save(result, band=1)
-
-        for idx, inlayer in enumerate(per_taxa):
-            assert inlayer.name
-            _, name = os.path.split(inlayer.name)
-            taxa = name[:-4]
-            species_count = int(total_counts[total_counts.taxa==taxa]["count"].values[0])
-            result._dataset.GetRasterBand(idx + 2).SetDescription(taxa)  # pylint: disable=W0212
-            scaled_filtered_layer = yo.where(
-                area_restore_filter != 0,
-                ((inlayer / area_restore_filter) * -1.0) / species_count,
-                float('nan')
-            )
-            scaled_filtered_layer.parallel_save(result, band=idx + 2)
+        bands.append(scaled_filtered_layer)
 
+    yg.to_geotiff(output_path, bands, labels, nodata=float('nan'))
 
 def main() -> None:
     parser = argparse.ArgumentParser(description="Scale final results for publication.")