tony2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################################\n",
    "# GeoHack 2023 Automated Fire Perimiter\n",
    "#\n",
    "#Provides functions to:\n",
    "# 1.1 > Classify Raster\n",
    "# 1.2 > Raster to Vector\n",
    "#     > ??? Project Vectors ???\n",
    "#     > Aggregate Polygons\n",
    "#     > Split\n",
    "#     > Combine\n",
    "#\n",
    "#Information:\n",
    "#     Automated Fire Perimeter model using classification and\n",
    "#     transformation techniques to rapidly build a fire\n",
    "#     shape file. This model reduces noise and accounts for\n",
    "#     layers of human subjectivity based on industry defined\n",
    "#     threshold. Prectical applications for fire inteligence\n",
    "#     are situational awareness of fire growth in a near \n",
    "#     real time environment. Each succesive update offers \n",
    "#     a POI for AI fire pretiction models from other platforms\n",
    "# \n",
    "# Acknowledgments: GeoHack 2023 Team 1 (\"Arno Roeder (Germany); Ricardo Lopes (Portugal); Sara Houser (Switzerland); \n",
    "#                  Flavia Ackermann (Switzerland); Curtis Doty (United States); ChatGPT (Open AI)\n",
    "#          > Info: Special thanks to the team for supporting contibutions which were instrumental to the build.\n",
    "#   \n",
    "# By: Tony Ramos\n",
    "# Date: 03/10/2023\n",
    "#\n",
    "############################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import rasterio\n",
    "import shapely\n",
    "import numpy as np\n",
    "import geopandas\n",
    "import scipy\n",
    "\n",
    "from matplotlib import pyplot\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# File Path\n",
    "\n",
    "work_dir = \"./data/\"\n",
    "image = os.path.join(work_dir, \"LWIR_QuickMosaic_16-bit_9327.tiff\")\n",
    "poly = os.path.join(work_dir, \"HeatPoly.shp\")\n",
    "# work_dir =\"Data\"\n",
    "# MillsFire=\"C:/Users/centr/GeoHack23-WildfirePerimeter/Data/LWIR_QuickMosaic_16-bit_9327.tiff\"\n",
    "OutputFilePath =\"./data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.1\n",
    "# CLASSIFY RASTER\n",
    "\n",
    "# Here's a step-by-step breakdown of the code:\n",
    "\n",
    "# The code defines a threshold value KneeThresh which is used to classify the raster data later on.\n",
    "\n",
    "# The with statement opens the input raster file using rasterio, and the raster data and metadata are \n",
    "# read and stored in raster_data and raster_meta variables, respectively.\n",
    "\n",
    "# The numpy library is used to classify the raster data based on the threshold value. The np.where() \n",
    "# function is used to create a new array where all values less than KneeThresh are set to 0, and all \n",
    "# other values are set to 1. The resulting array is stored in the classified_data variable.\n",
    "\n",
    "# The metadata of the input raster is updated to reflect the new data type and nodata value. \n",
    "# The rasterio.uint8 data type is used and nodata is set to None.\n",
    "\n",
    "# The output raster file is created using the rasterio.open() function. The output file is named \n",
    "# RasterClass.tif and is located in the OutputFilePath directory. The metadata of the input raster \n",
    "# is passed to the function as **raster_meta which unpacks the dictionary into keyword arguments. \n",
    "# The w mode is used to open the file for writing.\n",
    "\n",
    "# The dst.write() function writes the classified data to the output raster file. The astype() method \n",
    "# is used to convert the data type of the classified_data array to rasterio.uint16, which matches the \n",
    "# data type of the output raster file. The 1 parameter indicates that the data should be written to \n",
    "# the first band of the raster file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1.1\n",
    "# Open the raster file using rasterio\n",
    "KneeThresh = 33332\n",
    "with rasterio.open(image) as src:\n",
    "    # Read the raster data and metadata\n",
    "    raster_data = src.read(1)  # Read the first band of the raster\n",
    "    raster_meta = src.meta\n",
    "\n",
    "    # Classify the raster data\n",
    "    classified_data = np.where(raster_data < KneeThresh, 0, 1)\n",
    "\n",
    "# Update the metadata with the new data type and nodata value\n",
    "raster_meta.update(dtype=rasterio.uint8, nodata=None)\n",
    "\n",
    "# Write the classified raster to a new file\n",
    "with rasterio.open(OutputFilePath + 'RasterClass.tif', 'w', **raster_meta) as dst:\n",
    "    dst.write(classified_data.astype(rasterio.uint16), 1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.2\n",
    "# RASTER TO VECTOR\n",
    "\n",
    "# The with statement opens the input raster file using rasterio, and the raster data and metadata are read \n",
    "# and stored in raster_data and raster_meta variables, respectively.\n",
    "\n",
    "# The numpy library is used to classify the raster data based on a threshold value. The np.where() function \n",
    "# is used to create a new array where all values less than KneeThresh are set to 0, and all other values are set to 1. \n",
    "# The resulting array is stored in the classified_data variable.\n",
    "\n",
    "# The metadata of the input raster is updated to reflect the new data type and nodata value. The rasterio.int16 \n",
    "# data type is used and nodata is set to None.\n",
    "\n",
    "# The output raster file is created using the rasterio.open() function. The output file is named RasterClass.tif \n",
    "# and is located in the OutputFilePath directory. The metadata of the input raster is passed to the function as \n",
    "# **raster_meta which unpacks the dictionary into keyword arguments. The w mode is used to open the file for writing.\n",
    "\n",
    "# The dst.write() function writes the classified data to the output raster file. The 1 parameter indicates that the \n",
    "# data should be written to the first band of the raster file.\n",
    "\n",
    "# The features.shapes() function from the rasterio.features module is used to convert the classified raster to polygons. \n",
    "# The transform parameter is used to provide the affine transformation matrix that maps the raster coordinates to real-world \n",
    "# coordinates.\n",
    "\n",
    "# The fiona library is used to write the polygons to a shapefile. The output shapefile is named HeatPoly.shp \n",
    "# and is located in the OutputFilePath directory. The w mode is used to open the file for writing, and the \n",
    "# crs parameter is used to set the coordinate reference system (CRS) of the shapefile. The schema parameter is \n",
    "# used to define the schema of the shapefile, which in this case only includes a geometry field and no attribute fields.\n",
    "\n",
    "# A loop is used to iterate over each polygon in the shapes list. The loop checks if the polygon has a nonzero value, \n",
    "# indicating that it belongs to a class, and creates a new feature dictionary with the polygon geometry and no attribute \n",
    "# values. The dst.write() function writes the feature to the output shapefile.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1.2\n",
    "import rasterio\n",
    "from rasterio import features\n",
    "import fiona\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Open the raster file using rasterio\n",
    "with rasterio.open(image) as src:\n",
    "    # Read the raster data and metadata\n",
    "    raster_data = src.read(1)  # Read the first band of the raster\n",
    "    raster_meta = src.meta\n",
    "\n",
    "    # Classify the raster data\n",
    "    classified_data = np.where(raster_data < KneeThresh, 0, 1)\n",
    "\n",
    "# Update the metadata with the new data type and nodata value\n",
    "raster_meta.update(dtype=rasterio.int16, nodata=None)\n",
    "\n",
    "# Write the classified raster to a new file\n",
    "with rasterio.open(OutputFilePath + 'RasterClass.tif', 'w', **raster_meta) as dst:\n",
    "    dst.write(classified_data, 1)\n",
    "\n",
    "# Convert the classified raster to polygons\n",
    "shapes = features.shapes(classified_data, transform=raster_meta['transform'])\n",
    "\n",
    "# Write the polygons to a shapefile\n",
    "with fiona.open(OutputFilePath + 'HeatPoly.shp', 'w', 'ESRI Shapefile',crs=fiona.crs.from_epsg(4326), schema={'geometry': 'Polygon', 'properties': {}}) as dst:\n",
    "    for shape in shapes:\n",
    "        value = shape[1]\n",
    "        if value > 0:\n",
    "            feature = {'geometry': shape[0], 'properties': {}}\n",
    "            dst.write(feature)\n",
    "\n",
    "# # (Optional) Display the polygons using matplotlib\n",
    "# with fiona.open(OutputFilePath + 'HeatPoly.shp', 'r') as src:\n",
    "#     for i, layer in enumerate(src):\n",
    "#         if i == 0:\n",
    "#             x, y = [], []\n",
    "#             for poly in layer['geometry']['coordinates']:\n",
    "#                 x += [pt[0] for pt in poly]\n",
    "#                 y += [pt[1] for pt in poly]\n",
    "#             plt.plot(x, y, 'k.')\n",
    "#             plt.axis('equal')\n",
    "#             plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "\n",
    "# Define the source and target CRSs\n",
    "src_crs = 'EPSG:4326'  # EPSG:4326 is the default CRS for the input shapefile\n",
    "target_crs = 'EPSG:26910' # New projection NAD 1983 UTM Zone 10N\n",
    "\n",
    "# Read the input shapefile into a GeoDataFrame\n",
    "gdf = gpd.read_file(OutputFilePath + 'HeatPoly.shp', crs=src_crs)\n",
    "\n",
    "# Reproject the GeoDataFrame to the target CRS\n",
    "gdf = gdf.to_crs(target_crs)\n",
    "\n",
    "# Write the reprojected GeoDataFrame to a new shapefile\n",
    "gdf.to_file(OutputFilePath + 'HeatPoly_reprojected.shp', driver='ESRI Shapefile')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as shp\n",
    "from shapely.ops import unary_union\n",
    "from shapely.geometry import MultiPolygon\n",
    "\n",
    "# Load your polygons into a GeoDataFrame\n",
    "polygons = shp.read_file(work_dir + 'HeatPoly_reprojected.shp')\n",
    "\n",
    "# Set the distance within which you want to aggregate polygons\n",
    "distance = 40  # in meters\n",
    "\n",
    "# Create a buffer around each polygon using the distance\n",
    "buffered = polygons.geometry.buffer(distance)\n",
    "\n",
    "# Group buffered polygons that intersect with each other\n",
    "groups = buffered.unary_union\n",
    "\n",
    "# Convert the grouped polygons back to a GeoDataFrame\n",
    "# Convert the grouped polygons back to a GeoDataFrame\n",
    "if isinstance(groups, MultiPolygon):\n",
    "    polygons_list = [polygon for polygon in groups.geoms]\n",
    "else:\n",
    "    polygons_list = [groups]\n",
    "\n",
    "grouped_polygons = gpd.GeoDataFrame(\n",
    "    {'geometry': polygons_list},\n",
    "    crs=polygons.crs\n",
    ")\n",
    "\n",
    "# Save the aggregated polygons to a shapefile\n",
    "grouped_polygons.to_file(OutputFilePath + 'AggrigateEx1.shp')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the feature class\n",
    "fc = gpd.read_file(OutputFilePath + 'AggrigateEx1.shp')\n",
    "\n",
    "# Define a minimum area threshold for interior polygons (in square units of the CRS)\n",
    "threshold = 500\n",
    "\n",
    "# Iterate over each polygon in the feature class\n",
    "for index, row in fc.iterrows():\n",
    "    # Get the exterior and interior polygons\n",
    "    exterior = row.geometry.exterior\n",
    "    interiors = row.geometry.interiors\n",
    "    \n",
    "    # Filter out small interior polygons\n",
    "    interiors_filtered = [interior for interior in interiors if interior.area > threshold]\n",
    "    \n",
    "    # Create a new polygon with the filtered interiors\n",
    "    filtered_polygon = type(row.geometry)(exterior, interiors_filtered)\n",
    "    \n",
    "    # Replace the original geometry with the filtered geometry\n",
    "    fc.loc[index, 'geometry'] = filtered_polygon\n",
    "    \n",
    "# Save the updated feature class\n",
    "fc.to_file(OutputFilePath + 'Buff.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load your polygons into a GeoDataFrame\n",
    "polygons = gpd.read_file(work_dir + 'Buff.shp')\n",
    "\n",
    "# Set the distance within which you want to aggregate polygons\n",
    "distance = -43  # in meters\n",
    "\n",
    "# Create a buffer around each polygon using the distance\n",
    "buffered = polygons.geometry.buffer(distance)\n",
    "\n",
    "# Group buffered polygons that intersect with each other\n",
    "groups = buffered.unary_union\n",
    "\n",
    "# Convert the grouped polygons back to a GeoDataFrame\n",
    "# Convert the grouped polygons back to a GeoDataFrame\n",
    "if isinstance(groups, MultiPolygon):\n",
    "    polygons_list = [polygon for polygon in groups.geoms]\n",
    "else:\n",
    "    polygons_list = [groups]\n",
    "\n",
    "grouped_polygons = gpd.GeoDataFrame(\n",
    "    {'geometry': polygons_list},\n",
    "    crs=polygons.crs\n",
    ")\n",
    "\n",
    "# Save the aggregated polygons to a shapefile\n",
    "grouped_polygons.to_file(OutputFilePath + 'NegBuff.shp')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the polygon shapefile\n",
    "polygons = gpd.read_file(OutputFilePath + 'NegBuff.shp')\n",
    "\n",
    "# Dissolve the polygons based on a column called 'column_name'\n",
    "dissolved_polygons = polygons.dissolve(by='FID')\n",
    "\n",
    "# Save the dissolved polygons to a new shapefile\n",
    "dissolved_polygons.to_file(OutputFilePath + 'Dissolve.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Read in the feature class\n",
    "# fc = gpd.read_file(OutputFilePath + 'Union.shp')\n",
    "\n",
    "# # Define a minimum area threshold for interior polygons (in square units of the CRS)\n",
    "# threshold = 1\n",
    "\n",
    "# # Iterate over each polygon in the feature class\n",
    "# for index, row in fc.iterrows():\n",
    "#     # Get the exterior and interior polygons\n",
    "#     exterior = row.geometry.exterior\n",
    "#     interiors = row.geometry.interiors\n",
    "    \n",
    "#     # Filter out small interior polygons\n",
    "#     interiors_filtered = [interior for interior in interiors if interior.area > threshold]\n",
    "    \n",
    "#     # Create a new polygon with the filtered interiors\n",
    "#     filtered_polygon = type(row.geometry)(exterior, interiors_filtered)\n",
    "    \n",
    "#     # Replace the original geometry with the filtered geometry\n",
    "#     fc.loc[index, 'geometry'] = filtered_polygon\n",
    "    \n",
    "# # Save the updated feature class\n",
    "# fc.to_file(OutputFilePath + 'Union2.shp')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "wildfire",
   "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.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}