feat: generate h3 polygons and create webmap

.env.example

@@ -1 +1,2 @@
-MAPILLARY_CLIENT_TOKEN = "MY_MAPILLARY_CLIENT_TOKEN"
+MAPILLARY_CLIENT_TOKEN = "MY_MAPILLARY_CLIENT_TOKEN"
+MAPTILER_API_KEY = "MY_MAPTILER_API_KEY"

README.md

@@ -47,6 +47,7 @@ If you are interested in joining the project, please check out [`CONTRIBUTING.md
     - For example, download [`Three_Rivers_Michigan_USA_line.zip`](https://drive.google.com/file/d/1fpI4I5KP2WyVD5PeytW_hoXZswOt0dwA/view?usp=drive_link) to `data/raw/Three_Rivers_Michigan_USA_line.zip`. Note that this Google Drive link is only accessible to approved project members.
 4. Make a copy of the `.env.example` file, removing the `.example` from the end of the filename.
     - To download images from [Mapillary](https://www.mapillary.com/) you will need to create a (free) account and replace `MY_MAPILLARY_CLIENT_TOKEN` in the `.env` file with your own token. See the "Setting up API access and obtaining a client token" section on this [Mapillary help page](https://help.mapillary.com/hc/en-us/articles/360010234680-Accessing-imagery-and-data-through-the-Mapillary-API). You only need to enable READ access scope on your token.
+    - To use OpenStreetMap as a basemap for any webmaps generated, you will need a MapTiler API key. To get a free API key, follow the instructions on [this page](https://docs.maptiler.com/cloud/api/authentication-key/). Once you have an API key, add a new line to your `.env` file: `MAPTILER_API_KEY = "MY_MAPTILER_API_KEY"` - replace the text between the quotes with the key generated on the MapTiler account page. 
 
 ### 1. Sample points from roads data
 
@@ -104,6 +105,36 @@ saves an output to a new file.
 python -m src.assign_gvi_to_points data/raw/mapillary data/interim/Three_Rivers_Michigan_USA_points_images.gpkg data/processed/Three_Rivers_GVI.gpkg
 ```
 
+### 4. Visualize the results
+
+We provide an option here but we encourage you to explore different ways to visualize the results. We would love to know what you try that works and what doesn't. How can we improve on these guidance materials?
+
+#### Generate an H3 polygon layer
+
+We can generate an H3 polygon layer from the point layer. As an overlay it may make spatial trends more visible by merging some of the values from close together points. 
+
+### Example
+
+```bash
+# python -m src.create_webmap path/to/input_file.gpkg path/to/output_file.gpkg cell_resolution 
+python -m src.points_to_h3 data/processed/Three_Rivers_GVI.gpkg data/processed/Three_Rivers_h3_polygons_10.gpkg 10
+```
+
+The larger the number for the [H3 cell resolution](https://h3geo.org/docs/core-library/restable/), the smaller the individual hexagons. 
+
+#### Generate a web map
+
+We can generate an HTML file that contains javascript to display a web map showing the H3 polygons, styled by the mean GVI score for each polygon.
+
+To display an OpenStreetMap basemap under the data, you will need an API key from [MapTiler](https://www.maptiler.com/), a vector tiles provider. Once you have an API key, add it to your `.env` file (see [Setup section](#0-setup) for how to do this).
+
+### Example
+
+```bash
+# python -m src.create_webmap path/to/input_file.gpkg path/to/output/output_file.html default_zoom_for_webmap 
+python -m src.create_webmap data/processed/Three_Rivers_h3_polygons_10.gpkg data/processed/Three_Rivers_gvi_webmap.html 10
+```
+
 ## Config files
 
 > ![NOTE]

pyproject.toml

@@ -17,6 +17,7 @@ dependencies = [
   "folium",
   "geopandas",
   "geopy",
+  "h3pandas",
   "loguru",
   "mapclassify",
   "matplotlib",

src/create_webmap.py

@@ -0,0 +1,162 @@
+import os
+from pathlib import Path
+
+from dotenv import load_dotenv
+import geopandas as gpd
+from jinja2 import Environment, FileSystemLoader
+import matplotlib
+import numpy as np
+from shapely.geometry import box
+import typer
+
+try:
+    from typing import Annotated
+except ImportError:
+    # for Python 3.9
+    from typing_extensions import Annotated
+
+
+app = typer.Typer()
+
+
+def h3_to_webmap():
+    return
+
+
+@app.command()
+def main(
+    input_file: Annotated[
+        Path,
+        typer.Argument(help="Path to file containing h3 polygons."),
+    ],
+    filename: Annotated[
+        str, typer.Argument(help="(Optional) Path to file for HTML output.")
+    ] = "./data/processed/gvi_webmap.html",
+    zoom_level: Annotated[
+        int,
+        typer.Argument(
+            help="""(Optional) Starting zoom level for webmap. 
+			Takes integer between 0 (small scale) and 20 (large scale).
+			"""
+        ),
+    ] = 12,
+):
+    gdf = gpd.read_file(input_file)
+    # if crs is not 4326, convert to 4326
+    if gdf.crs == "EPSG:4326":
+        pass
+    else:
+        gdf = gdf.to_crs("EPSG:4326")
+
+    # Round GVI score to make map labels more readable
+    gdf["gvi_score"] = round(gdf["gvi_score"], 2)
+
+    # get central coordinates of all features
+    centre = (
+        str(gdf.dissolve().centroid.x.values[0])
+        + ", "
+        + str(gdf.to_crs("4326").dissolve().centroid.y.values[0])
+    )
+    centre_str = "[" + centre + "]"
+
+    # Calculate datasdet bounds (with a buffer of 0.5 degrees) to limit webmap bounds
+    gdf_bounds = gdf.total_bounds
+    bbox = box(*gdf_bounds).buffer(0.25, cap_style="square", join_style="mitre")
+    bounds = bbox.bounds
+    bounds_str = (
+        "["
+        + str(bounds[0])
+        + ","
+        + str(bounds[1])
+        + "], ["
+        + str(bounds[2])
+        + ","
+        + str(bounds[3])
+        + "]"
+    )
+
+    # Load API key for map tiles from environment variable
+    load_dotenv()
+
+    if "MAPTILER_API_KEY" in os.environ:
+        maptiler_api_key = os.getenv("MAPTILER_API_KEY")
+    else:
+        maptiler_api_key = "None"
+
+    if maptiler_api_key in ["MY_MAPTILER_API_KEY", ""]:
+        maptiler_api_key = "None"
+
+    # Lookup the colourmap values for each GVI score
+    cmap = matplotlib.colormaps["Greens"]
+    gdf["gvi_norm"] = (gdf.gvi_score - np.min(gdf.gvi_score)) / (
+        np.max(gdf.gvi_score) - np.min(gdf.gvi_score)
+    )
+    gdf["html_color"] = gdf["gvi_norm"].apply(
+        lambda x: matplotlib.colors.rgb2hex(cmap(x))
+    )
+
+    # Generate divs for legend
+    # Pick 4 evenly-spaced values from the gvi scores to use in the legend
+    legend_gvi = list(
+        np.arange(
+            gdf.gvi_score.min(),
+            gdf.gvi_score.max(),
+            (gdf.gvi_score.max() - gdf.gvi_score.min()) / 4,
+            dtype=int,
+        )
+    )
+
+    # Generate labels by looking up what the GVI score would be for those values
+    legend_label_1 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[0]
+    )
+    legend_label_2 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[33]
+    )
+    legend_label_3 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[66]
+    )
+    legend_label_4 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[99]
+    )
+
+    # Normalise the label values to lookup against the colourmap
+    legend_gvi_norm = (legend_gvi - np.min(legend_gvi)) / (
+        np.max(legend_gvi) - np.min(legend_gvi)
+    )
+
+    # Generate the html colour code from the normalised values
+    legend_colours = []
+    for i in legend_gvi_norm:
+        legend_colours.append(matplotlib.colors.rgb2hex(cmap(i)))
+    # Assign patch colours to use in HTML template
+    legend_patch_1, legend_patch_2, legend_patch_3, legend_patch_4 = legend_colours
+
+    # Load the MapLibre HMTL template
+    environment = Environment(loader=FileSystemLoader("src/templates"))
+    template = environment.get_template("maplibre_template.html")
+
+    # Generate the HTML file from the template, filling dynamic values
+    with open(filename, mode="w", encoding="utf-8") as message:
+        message.write(
+            template.render(
+                title="GVI score hex map",
+                geojson=gdf.to_json(),
+                centre_coords=centre_str,
+                zoom=zoom_level,
+                bounds=bounds_str,
+                maptiler_api_key=maptiler_api_key,
+                legend_label_1=legend_label_1,
+                legend_label_2=legend_label_2,
+                legend_label_3=legend_label_3,
+                legend_label_4=legend_label_4,
+                legend_patch_1=legend_patch_1,
+                legend_patch_2=legend_patch_2,
+                legend_patch_3=legend_patch_3,
+                legend_patch_4=legend_patch_4,
+            )
+        )
+
+
+if __name__ == "__main__":
+    app()

src/points_to_h3.py

@@ -0,0 +1,94 @@
+import os
+from pathlib import Path
+
+import geopandas as gpd
+import h3pandas
+import typer
+
+try:
+    from typing import Annotated
+except ImportError:
+    # for Python 3.9
+    from typing_extensions import Annotated
+
+app = typer.Typer()
+
+
+@app.command()
+def main(
+    input_file: Annotated[
+        Path,
+        typer.Argument(help="Path to file containing point layer with GVI scores."),
+    ],
+    output_file: Annotated[
+        Path,
+        typer.Argument(
+            help="File to write output data to (can specify any GDAL-supported format)"
+        ),
+    ],
+    cell_resolution: Annotated[
+        int,
+        typer.Argument(
+            help="""H3 cell resolution to aggregate to, 
+            between 0 (largest) and 15 (smallest)"""
+        ),
+    ] = 10,
+):
+    """
+    Aggregates points to h3 hex cells.
+
+    Args:
+            input_file: Path to file containing point layer with GVI scores.
+            cell_resolution: H3 cell resolution to aggregate to,
+                between 0 (largest) and 15 (smallest)
+            aggregation_operations:
+            output_file: File to write output data to
+                (can specify any GDAL-supported format)
+
+    Returns:
+            File containing h3 polygons with aggregated GVI scores
+
+    """
+    # Check input file exists
+    if os.path.exists(input_file):
+        pass
+    else:
+        raise ValueError("Input file could not be found")
+
+    # Check input file is a valid file for GeoPandas
+    try:
+        gpd.read_file(input_file)
+    except Exception as e:
+        raise e
+
+    # Check data contains point features
+    if "Point" in gpd.read_file(input_file).geometry.type.unique():
+        pass
+    else:
+        raise Exception("Expected point data in interim data file but none found")
+
+    # Check data contains numeric gvi_score field
+
+    # Load input data
+    gdf = gpd.read_file(input_file)
+
+    # Exclude points with no GVI score
+    gdf = gdf[~gdf.gvi_score.isna()]
+
+    # Assign points to h3 cells at the selected resolution
+    gdf_h3 = gdf.h3.geo_to_h3(cell_resolution).reset_index()
+
+    # Aggregate the points to the assigned h3 cell
+    gvi_mean = gdf_h3.groupby("h3_" + f"{cell_resolution:02}").agg(
+        {"gvi_score": "mean"}
+    )
+
+    # Convert the h3 cells to polygons
+    gvi_hex = gvi_mean.h3.h3_to_geo_boundary()
+
+    # Export the h3 polygons to the specified output file
+    gvi_hex.to_file(output_file)
+
+
+if __name__ == "__main__":
+    app()