TANABBO_forest_edge.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.merge import merge\n",
    "from scipy.ndimage import label, binary_erosion\n",
    "from shapely.geometry import Point, MultiPoint\n",
    "from shapely.ops import nearest_points, transform\n",
    "import geopandas as gpd\n",
    "import pyproj\n",
    "import glob\n",
    "from rasterio import features\n",
    "from geopandas.tools import sjoin\n",
    "from shapely.geometry import shape\n",
    "from scipy.ndimage import binary_dilation\n",
    "from shapely.ops import transform as shapely_transform\n",
    "from rasterio.transform import from_origin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge yearly forest loss layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_raster_layers(file_pattern, output_file):\n",
    "    file_list = glob.glob(file_pattern)\n",
    "    mosaic_list = []\n",
    "    out_trans = None\n",
    "\n",
    "    for fp in file_list:\n",
    "        with rasterio.open(fp) as src:\n",
    "            mosaic, transform = merge([src])\n",
    "            mosaic_list.append(mosaic)\n",
    "            if out_trans is None:\n",
    "                out_trans = transform\n",
    "\n",
    "    # Combine the mosaics into a single array\n",
    "    final_mosaic = np.concatenate(mosaic_list, axis=0)\n",
    "\n",
    "    # Convert to binary\n",
    "    final_mosaic_binary = (final_mosaic > 0).astype(rasterio.uint8)\n",
    "\n",
    "    # Assuming all bands are identical, just keep the first band\n",
    "    final_mosaic_binary = final_mosaic_binary[0:1, :, :]\n",
    "\n",
    "    out_meta = src.meta.copy()\n",
    "    out_meta.update({\n",
    "        \"height\": final_mosaic_binary.shape[1],\n",
    "        \"width\": final_mosaic_binary.shape[2],\n",
    "        \"transform\": out_trans,\n",
    "        \"count\": 1,  # Set band count to 1\n",
    "        \"dtype\": 'uint8',  # Set data type to uint8 for binary data\n",
    "        \"nodata\": 0  # Set nodata value suitable for uint8 type\n",
    "    })\n",
    "\n",
    "    with rasterio.open(output_file, \"w\", **out_meta) as dest:\n",
    "        dest.write(final_mosaic_binary)\n",
    "\n",
    "# Example of usage\n",
    "merge_raster_layers(\"bb_spot_20[1][1-8].tif\", \"bb_spot_merged.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Update forest mask "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_raster_data(file_path):\n",
    "    with rasterio.open(file_path) as src:\n",
    "        return src.read(1), src.transform, src.crs\n",
    "\n",
    "bb_spot_merged, forest_transform, forest_crs = load_raster_data(\"bb_spot_merged.tif\")\n",
    "s50mask, _, _ = load_raster_data(\"s50mask.tif\")\n",
    "\n",
    "def update_mask(merged_raster_file, mask_file, output_file):\n",
    "    with rasterio.open(merged_raster_file) as src:\n",
    "        merged_raster = src.read(1)  # Read the first band of bb_spot_merged.tif\n",
    "        crs = src.crs\n",
    "        transform = src.transform\n",
    "\n",
    "    with rasterio.open(mask_file) as src:\n",
    "        mask = src.read(1)  # Read the first band of s50mask.tif\n",
    "\n",
    "    # Wherever merged_raster is 1, set mask to 0\n",
    "    updated_mask = np.where(merged_raster == 1, 0, mask)\n",
    "\n",
    "    # Update metadata for output file\n",
    "    out_meta = src.meta.copy()\n",
    "    out_meta.update({\n",
    "        \"dtype\": str(updated_mask.dtype),\n",
    "        \"nodata\": 0  # Use 0 as the nodata value\n",
    "    })\n",
    "\n",
    "    with rasterio.open(output_file, \"w\", **out_meta) as dest:\n",
    "        dest.write(updated_mask, 1)\n",
    "\n",
    "# Update mask\n",
    "update_mask(\"bb_spot_merged.tif\", \"s50mask.tif\", \"s50mask_upd.tif\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forest edge mask generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def identify_and_filter_edges(bb_spot_merged, s50mask):\n",
    "    labeled_array, num_features = label(bb_spot_merged)\n",
    "\n",
    "    edge = np.zeros_like(labeled_array, dtype=np.bool_)\n",
    "    for i in range(1, num_features + 1):\n",
    "        object_mask = (labeled_array == i)\n",
    "        object_border = object_mask & np.logical_not(binary_erosion(object_mask))\n",
    "        edge = np.logical_or(edge, object_border)\n",
    "\n",
    "    edge = edge & bb_spot_merged.astype(np.bool_) & (bb_spot_merged > 0)\n",
    "    labeled_edge, num_edge_features = label(edge)\n",
    "\n",
    "    forest_edge_mask = np.zeros_like(labeled_edge, dtype=np.bool_)\n",
    "    for i in range(1, num_edge_features + 1):\n",
    "        object_mask = (labeled_edge == i).astype(bool)\n",
    "        if np.any(object_mask & s50mask.astype(bool)):\n",
    "            forest_edge_mask = np.logical_or(forest_edge_mask, object_mask)\n",
    "\n",
    "    # Using hardcoded transform and crs values from the loaded raster\n",
    "    hardcoded_transform = forest_transform\n",
    "    hardcoded_crs = forest_crs\n",
    "\n",
    "    # Save the forest_edge_mask as a GeoTIFF using the hardcoded values\n",
    "    with rasterio.open(\"forest_edge_mask.tif\", 'w', driver='GTiff', height=forest_edge_mask.shape[0],\n",
    "                       width=forest_edge_mask.shape[1], count=1, dtype=rasterio.uint8, crs=hardcoded_crs, transform=hardcoded_transform) as dst:\n",
    "        dst.write((forest_edge_mask * 255).astype(rasterio.uint8), 1)\n",
    "\n",
    "    return forest_edge_mask, edge\n",
    "\n",
    "forest_edge_mask, edge = identify_and_filter_edges(bb_spot_merged, s50mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pheromone trap installation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\trubin\\.pyenv\\pyenv-win\\versions\\3.11.1\\Lib\\site-packages\\IPython\\core\\interactiveshell.py:3505: FutureWarning: The `op` parameter is deprecated and will be removed in a future release. Please use the `predicate` parameter instead.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    }
   ],
   "source": [
    "def generate_and_filter_points(forest_edge_mask, s50mask, forest_transform, forest_crs):\n",
    "    rows, cols = np.where(forest_edge_mask == 1)\n",
    "    geometries_mask = [Point(forest_transform * (col + 0.5, row + 0.5)) for col, row in zip(cols, rows)]\n",
    "    gdf_mask = gpd.GeoDataFrame(geometry=geometries_mask)\n",
    "    gdf_mask.set_crs(forest_crs, inplace=True)\n",
    "\n",
    "    mask_polygon_features = [shape(shape_val) for shape_val, val in features.shapes(s50mask.astype(np.uint8), transform=forest_transform) if val == 1]\n",
    "\n",
    "    if not mask_polygon_features:\n",
    "        raise ValueError(\"No valid geometries found in s50mask!\")\n",
    "\n",
    "    mask_gdf = gpd.GeoDataFrame({'geometry': mask_polygon_features})\n",
    "    mask_gdf.crs = forest_crs\n",
    "\n",
    "    gdf_mask_within = sjoin(gdf_mask, mask_gdf, how=\"inner\", op=\"within\")\n",
    "    all_centroids_mask_filtered = gdf_mask_within['geometry'].values\n",
    "    \n",
    "    return all_centroids_mask_filtered\n",
    "\n",
    "# Generate and filter points\n",
    "all_centroids_mask_filtered = generate_and_filter_points(forest_edge_mask, s50mask, forest_transform, forest_crs)\n",
    "\n",
    "# Filter the centroids based on the 50m distance criterion\n",
    "def determine_utm_zone(centroid):\n",
    "    zone_number = (int((centroid.x + 180) / 6) % 60) + 1\n",
    "    if centroid.y > 0:\n",
    "        return f\"EPSG:326{zone_number:02}\"  # Northern hemisphere\n",
    "    else:\n",
    "        return f\"EPSG:327{zone_number:02}\"  # Southern hemisphere\n",
    "\n",
    "def generate_filtered_centroids_shp(edge, forest_transform, forest_crs, all_centroids_mask_filtered, output_file_path):\n",
    "    \"\"\"\n",
    "    Generate a shapefile of filtered centroids based on the provided edge and mask raster data.\n",
    "\n",
    "    Parameters:\n",
    "        edge (np.array): A 2D numpy array representing the edge raster.\n",
    "        forest_transform (Affine): Affine transformation for converting pixel to spatial coordinates.\n",
    "        forest_crs (CRS): Coordinate reference system for spatial data.\n",
    "        all_centroids_mask_filtered (list of shapely.Point): List of filtered centroid points.\n",
    "        output_file_path (str): File path to save the resulting shapefile.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Generate all_centroids.shp based on the edge\n",
    "    rows, cols = np.where(edge == 1)\n",
    "    geometries = [Point(forest_transform * (col + 0.5, row + 0.5)) for col, row in zip(cols, rows)]\n",
    "    gdf = gpd.GeoDataFrame(geometry=geometries)\n",
    "    gdf.set_crs(forest_crs, inplace=True)\n",
    "\n",
    "    # List to hold all centroid coordinates\n",
    "    all_centroids = [(geometry.x, geometry.y) for geometry in gdf.geometry]\n",
    "\n",
    "    # Get UTM Zone\n",
    "    utm_zone = determine_utm_zone(Point(all_centroids[0]))\n",
    "\n",
    "    # Filter the centroids based on the 50m distance criterion\n",
    "    project_to_utm = pyproj.Transformer.from_crs(forest_crs, utm_zone, always_xy=True).transform\n",
    "    project_to_wgs84 = pyproj.Transformer.from_crs(utm_zone, forest_crs, always_xy=True).transform\n",
    "    all_centroids_mask_utm = [transform(project_to_utm, Point(c)) for c in all_centroids_mask_filtered]\n",
    "\n",
    "    # Filter based on distance criterion\n",
    "    final_centroids_mask_utm = []\n",
    "    for point in all_centroids_mask_utm:\n",
    "        nearest = nearest_points(point, MultiPoint(final_centroids_mask_utm)) if final_centroids_mask_utm else [None, Point(9999, 9999)]\n",
    "        if point.distance(nearest[1]) > 50:\n",
    "            final_centroids_mask_utm.append(point)\n",
    "    final_centroids_mask = [transform(project_to_wgs84, point) for point in final_centroids_mask_utm]\n",
    "\n",
    "    # Perform dilation on the final_centroids_mask\n",
    "    mask_shape = edge.shape\n",
    "    final_mask = np.zeros(mask_shape, dtype=bool)\n",
    "\n",
    "    for point in final_centroids_mask_utm:\n",
    "        col, row = ~forest_transform * point.coords[0]\n",
    "        if 0 <= row < mask_shape[0] and 0 <= col < mask_shape[1]:\n",
    "            final_mask[int(row), int(col)] = True\n",
    "\n",
    "    # Apply binary dilation and identify isolated regions\n",
    "    dilated_mask = binary_dilation(final_mask, structure=np.ones((3,3)))\n",
    "    \n",
    "    # [NOTE]: Replace `s50mask` with your actual mask array here. \n",
    "    isolated_mask = dilated_mask & np.logical_not(s50mask)\n",
    "    isolated_rows, isolated_cols = np.where(isolated_mask)\n",
    "    isolated_points = [forest_transform * (col + 0.5, row + 0.5) for col, row in zip(isolated_cols, isolated_rows)]\n",
    "    isolated_geoms = [Point(p) for p in isolated_points]\n",
    "\n",
    "    # Drop isolated points from final_centroids_mask\n",
    "    final_centroids_mask = [point for point in final_centroids_mask if point not in isolated_geoms]\n",
    "\n",
    "    # Save the filtered centroids to a shapefile\n",
    "    geometry_mask = [Point(point.x, point.y) for point in final_centroids_mask]\n",
    "    gdf_filtered_mask = gpd.GeoDataFrame(geometry=geometry_mask)\n",
    "    gdf_filtered_mask.set_crs(forest_crs, inplace=True)\n",
    "    #gdf_filtered_mask.to_file(output_file_path)\n",
    "\n",
    "generate_filtered_centroids_shp(edge, forest_transform, forest_crs, all_centroids_mask_filtered, \"filtered_centroids_mask_upd.shp\")\n",
    "\n",
    "def remove_single_points(s50mask_upd_filepath, bb_spot_merged_filepath, gdf_filtered_mask, output_filepath):\n",
    "    \"\"\"\n",
    "    Removes single points based on conditions applied to two raster masks and saves the \n",
    "    filtered centroids to a shapefile.\n",
    "\n",
    "    Parameters:\n",
    "        s50mask_upd_filepath (str): Filepath to the \"s50mask_upd.tif\" raster.\n",
    "        bb_spot_merged_filepath (str): Filepath to the \"bb_spot_merged.tif\" raster.\n",
    "        gdf_filtered_mask (GeoDataFrame): GeoDataFrame containing points to be potentially removed.\n",
    "        output_filepath (str): Filepath to store the resulting shapefile with removed points.\n",
    "    \"\"\"\n",
    "\n",
    "    with rasterio.open(s50mask_upd_filepath) as src:\n",
    "        s50mask_upd = src.read(1)\n",
    "        raster_transform = src.transform\n",
    "\n",
    "    with rasterio.open(bb_spot_merged_filepath) as src_bb:\n",
    "        bb_mask = src_bb.read(1)\n",
    "\n",
    "    # Function to check if a point's corresponding pixel and its surroundings meet the conditions\n",
    "    def should_drop(row, col):\n",
    "        condition1 = bb_mask[row][col] == 1 and all(\n",
    "            bb_mask[i][j] == 0 for i in range(row-1, row+2) for j in range(col-1, col+2) if (i, j) != (row, col)\n",
    "        )\n",
    "\n",
    "        condition2 = s50mask_upd[row][col] == 0 and all(\n",
    "            s50mask_upd[i][j] == 1 for i in range(row-1, row+2) for j in range(col-1, col+2) if (i, j) != (row, col)\n",
    "        )\n",
    "\n",
    "        return condition1 and condition2\n",
    "\n",
    "    # Find the centroids that meet the conditions and remove them\n",
    "    to_drop = []\n",
    "    for idx, point in gdf_filtered_mask.iterrows():\n",
    "        col, row = ~raster_transform * (point.geometry.x, point.geometry.y)\n",
    "        col, row = int(col), int(row)\n",
    "        if should_drop(row, col):\n",
    "            to_drop.append(idx)\n",
    "    gdf_filtered_mask = gdf_filtered_mask.drop(index=to_drop)\n",
    "    #gdf_filtered_mask.to_file(output_filepath)\n",
    "\n",
    "gdf_filtered_mask = gpd.read_file(\"filtered_centroids_mask_upd.shp\")\n",
    "remove_single_points(\"s50mask_upd.tif\", \"bb_spot_merged.tif\", gdf_filtered_mask, \"filtered_centroids_mask_upd_cleaned.shp\")\n",
    "\n",
    "\n",
    "# South feature\n",
    "def has_south_neighbor_of_value(row, col, raster, value):\n",
    "    \"\"\"\n",
    "    Check if the pixel south to (row, col) in 'raster' has the given 'value'.\n",
    "    \"\"\"\n",
    "    # Check if the south pixel is within the raster bounds\n",
    "    if row + 1 < raster.shape[0]:\n",
    "        return raster[row + 1, col] == value\n",
    "    return False\n",
    "\n",
    "def filter_centroids_without_south_neighbor(raster_path, input_shp_path, output_shp_path, value):\n",
    "    \"\"\"\n",
    "    Filters out points from a shapefile that have a southern neighbor with a specific value in a raster.\n",
    "\n",
    "    Parameters:\n",
    "        raster_path (str): Path to the raster file.\n",
    "        input_shp_path (str): Path to the input shapefile.\n",
    "        output_shp_path (str): Path to save the filtered shapefile.\n",
    "        value (int or float): The value to check in the southern neighbor in the raster.\n",
    "    \"\"\"\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        raster_data = src.read(1)\n",
    "        raster_transform = src.transform\n",
    "    \n",
    "    gdf = gpd.read_file(input_shp_path)\n",
    "    \n",
    "    to_keep = []\n",
    "    for _, point in gdf.iterrows():\n",
    "        col, row = ~raster_transform * (point.geometry.x, point.geometry.y)\n",
    "        col, row = int(col), int(row)\n",
    "        \n",
    "        if not has_south_neighbor_of_value(row, col, raster_data, value):\n",
    "            to_keep.append(True)\n",
    "        else:\n",
    "            to_keep.append(False)\n",
    "    \n",
    "    gdf_filtered = gdf.loc[to_keep]\n",
    "    gdf_filtered.to_file(output_shp_path)\n",
    "\n",
    "filter_centroids_without_south_neighbor(\n",
    "    \"s50mask_upd.tif\", \n",
    "    \"filtered_centroids_mask_upd_cleaned.shp\", \n",
    "    \"pheromone_trap.shp\", \n",
    "    1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Antiattractant installation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Buffer creation around the forest loss adjacent areas in the updated forest mask is completed.\n",
      "Filtering of forest_edge_buffer complete, resulting in antiattractant.tif.\n",
      "Number of pixels with value 1 in the final raster: 641\n"
     ]
    }
   ],
   "source": [
    "# Function to read a raster file\n",
    "def read_raster(file_path):\n",
    "    with rasterio.open(file_path) as src:\n",
    "        return src.read(1), src.profile\n",
    "\n",
    "# Function to write a raster file\n",
    "def write_raster(output_path, data, profile):\n",
    "    profile.update(\n",
    "        dtype=rasterio.uint8,\n",
    "        count=1,\n",
    "        compress='lzw')\n",
    "    with rasterio.open(output_path, 'w', **profile) as dst:\n",
    "        dst.write(data.astype(rasterio.uint8), 1)\n",
    "\n",
    "# Read the updated forest mask and forest loss rasters\n",
    "forest_mask, forest_profile = read_raster(\"s50mask_upd.tif\")\n",
    "loss_mask, _ = read_raster(\"bb_spot_merged.tif\")\n",
    "\n",
    "# Ensure the masks are boolean\n",
    "forest_mask = forest_mask.astype(bool)\n",
    "loss_mask = loss_mask.astype(bool)\n",
    "\n",
    "# Create a dilation of the loss mask\n",
    "dilated_loss_mask = binary_dilation(loss_mask, iterations=1)\n",
    "\n",
    "# Find the edge pixels: pixels that are in the dilated loss mask AND in the updated forest mask\n",
    "edge_pixels = dilated_loss_mask & forest_mask\n",
    "\n",
    "# Convert the edge pixels mask back to uint8 for raster output\n",
    "edge_pixels = edge_pixels.astype(np.uint8)\n",
    "\n",
    "# Write the result\n",
    "write_raster(\"forest_edge_buffer.tif\", edge_pixels, forest_profile)\n",
    "\n",
    "print(\"Buffer creation around the forest loss adjacent areas in the updated forest mask is completed.\")\n",
    "\n",
    "# Read the forest_edge_buffer raster\n",
    "forest_edge_buffer, buffer_profile = read_raster(\"forest_edge_buffer.tif\")\n",
    "\n",
    "# Read and merge all bands of bb_spot_merged.tif into a binary mask\n",
    "with rasterio.open(\"bb_spot_merged.tif\") as src:\n",
    "    # Read all bands and stack them into a single ndarray (assumes all bands are of the same dimensions)\n",
    "    bb_spot_merged = src.read().astype(bool)\n",
    "    # Reduce to a binary mask where if any band has a pixel value of True, it results in True\n",
    "    bb_spot_merged_combined = np.any(bb_spot_merged, axis=0)\n",
    "\n",
    "# Make sure we have a boolean array for the forest edge buffer\n",
    "forest_edge_buffer_bool = forest_edge_buffer.astype(bool)\n",
    "\n",
    "# Initialize the output mask with the current forest edge buffer\n",
    "antiattractant_mask = np.copy(forest_edge_buffer_bool)\n",
    "\n",
    "# Iterate over each pixel, skipping the first row to avoid index errors\n",
    "for row in range(1, forest_edge_buffer.shape[0]):\n",
    "    for col in range(forest_edge_buffer.shape[1]):\n",
    "        # If the pixel in forest_edge_buffer is True and the pixel to the north in bb_spot_merged_combined is also True\n",
    "        if forest_edge_buffer_bool[row, col] and bb_spot_merged_combined[row-1, col]:\n",
    "            antiattractant_mask[row, col] = False\n",
    "\n",
    "# Convert the boolean mask back to an integer type suitable for writing to a TIFF\n",
    "antiattractant = antiattractant_mask.astype(np.uint8)\n",
    "\n",
    "# Write the antiattractant mask to a new file\n",
    "write_raster(\"antiattractant.tif\", antiattractant, buffer_profile)\n",
    "\n",
    "print(\"Filtering of forest_edge_buffer complete, resulting in antiattractant.tif.\")\n",
    "\n",
    "# Counting pixels with value of 1\n",
    "num_pixels = np.sum(antiattractant == 1)\n",
    "print(f\"Number of pixels with value 1 in the final raster: {num_pixels}\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}