feat: generate h3 polygons and create webmap (#46)

* Create points_to_h3.py

* Revert "Create points_to_h3.py"

This reverts commit 9b33217.

* Create points_to_h3.py

* Updates to points_to_h3

- add h3pandas to dependencies
- add app launch to end of file
- use single aggregation (mean) instead of multiple
- change order of arguments based on defaults

* Add webmap code and template

* Tidying and formatting webmap script

- adding Typer help and defaults
- adding code to generate colours and breakpoints for map
- formatting

* Fix interpolation breaks

- fix interpolation breaks
- fix hard coded template path

* Fix code to generate polygon colours

Hard code colours into dataframe/json rather than using MapLibre interpolation

* Delete gvi_webmap.html

remove temp file

* Applying ruff checks

* Update maplibre_template.html

- add pop to show gvi score when polygon is clicked

* Update create_webmap.py

- Round GVI scores to 2 decimal places for neater pop-ups

* Generate legend from GVI scores

* Update create_webmap.py

- Reduce bounding box buffer from 0.5 degrees to 0.25
- Set default zoom to 12

* Update .env.example

Update sample .env file to include placeholder for Maptiler api key

* Update create_webmap.py

Fix line lengths

* Update points_to_h3.py

Formatting

* visualize step described in readme

* remove reference to folder not in repo

* alphabetize new dependency

* Update README.md

add info on MapTiler API key

* Fix single digit resolution error

Add a leading `0` to cell resolution if it is  a single digit

* Use empy basemap if no Maptiler key provided

The Jinja template will use an empty basemap if there is no Maptiler API key in the .env file

* Formatting

* fix env variable and update template

* Combine output path and filename

* docs: filename and path is one argument

* add logiv for missing env variable

* fix for blank env variable

* fix: ruff format

* fix: ruff fix import block is un-sorted or un-formatted

* fix: change web map colormap

* fix legend by changing int to float

* Ignore linter for import; needed to register accessor

---------

Co-authored-by: Dan Joseph <[email protected]>
Co-authored-by: Dan Joseph <[email protected]>
Co-authored-by: Alexei <[email protected]>
Co-authored-by: Jay Qi <[email protected]>

.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=float,
+        )
+    )
+
+    # 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], 1
+    )
+    legend_label_2 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[33], 1
+    )
+    legend_label_3 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[66], 1
+    )
+    legend_label_4 = round(
+        np.linspace(gdf.gvi_score.min(), gdf.gvi_score.max(), 100)[99], 1
+    )
+
+    # 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  # noqa: F401
+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()