jupyter: Raster support for interactive maps with folium (#1769)

* Add temporary WGS84-Pseudo-Mercator location for rasters
* Add `add_raster` method for `InteractiveMap` class
* Added resolution estimator for rasters (for speed, reproject only as many cells as needed)
* Added `save()` method for saving folium map as HTML file
* Change `GrassRenderer` PNG background to transparent
* Added `utils.py' for clarity in the interactive class
* Added rcfile cleanup
* Added temporary directory with `tempfile`
* Update demonstration notebook to include raster example

doc/notebooks/jupyter_integration.ipynb

@@ -51,16 +51,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Let's check our import by printing docstring for init()\n",
-    "gj.init?"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
     "# Start GRASS Session\n",
     "gj.init(\"../../data/grassdata\", \"nc_basic_spm_grass7\", \"user1\")"
    ]
@@ -147,24 +137,6 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### Error Handling"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# If we pass a non-display related module to GrassRenderer, it returns an error\n",
-    "\n",
-    "# r.run(\"r.watershed\", map=\"elevation\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
     "## Interactive Map Display"
    ]
   },
@@ -174,12 +146,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Demonstration of GrassRenderer for interactive map display with folium\n",
-    "\n",
+    "# Create Interactive Map\n",
     "m = gj.InteractiveMap(width = 600)\n",
-    "m.add_vector(\"boundary_region@PERMANENT\")\n",
+    "\n",
+    "# Add raster, vector and layer control to map\n",
+    "m.add_raster(\"elevation\")\n",
+    "m.add_vector(\"roadsmajor\")\n",
     "m.add_layer_control(position = \"bottomright\")\n",
-    "m.add_vector(\"schools\")\n",
+    "\n",
+    "# Display map\n",
     "m.show()"
    ]
   },

python/grass/jupyter/Makefile

@@ -8,7 +8,8 @@ DSTDIR = $(ETC)/python/grass/jupyter
 MODULES = \
 	setup \
 	display \
-	interact_display
+	interact_display \
+	utils
 
 PYFILES := $(patsubst %,$(DSTDIR)/%.py,$(MODULES) __init__)
 PYCFILES := $(patsubst %,$(DSTDIR)/%.pyc,$(MODULES) __init__)

python/grass/jupyter/__init__.py

@@ -1,3 +1,4 @@
 from .setup import *
 from .interact_display import *
 from .display import *
+from .utils import *

python/grass/jupyter/display.py

@@ -71,6 +71,7 @@ class GrassRenderer:
         self._env["GRASS_RENDER_TEXT_SIZE"] = str(text_size)
         self._env["GRASS_RENDER_IMMEDIATE"] = renderer
         self._env["GRASS_RENDER_FILE_READ"] = "TRUE"
+        self._env["GRASS_RENDER_TRANSPARENT"] = "TRUE"
 
         # Create PNG file for map
         # If not user-supplied, we will write it to a map.png in a

python/grass/jupyter/interact_display.py

@@ -6,60 +6,69 @@
 #
 # COPYRIGHT: (C) 2021 Caitlin Haedrich, and by the GRASS Development Team
 #
-#            This program is free software under the GNU Gernal Public
-#            License (>=v2). Read teh file COPYING that comes with GRASS
+#            This program is free software under the GNU General Public
+#            License (>=v2). Read the file COPYING that comes with GRASS
 #            for details.
 
 import os
+import sys
+import tempfile
+import weakref
 from pathlib import Path
-
 import folium
-
 import grass.script as gs
+from .display import GrassRenderer
+from .utils import (
+    estimate_resolution,
+    get_location_proj_string,
+    get_region,
+    reproject_region,
+    setup_location,
+)
 
 
 class InteractiveMap:
-    """This class creates interative GRASS maps with folium"""
+    """This class creates interative GRASS maps with folium.
+
+    Basic Usage:
+    >>> m = InteractiveMap()
+    >>> m.add_vector("streams")
+    >>> m.add_raster("elevation")
+    >>> m.add_layer_control()
+    >>> m.show()
+    """
 
     def __init__(self, width=400, height=400):
-        """This initiates a folium map centered on g.region.
+        """Creates a blank folium map centered on g.region.
 
-        Keyword arguments:
-            height -- height in pixels of figure (default 400)
-            width -- width in pixels of figure (default 400)"""
+        :param int height: height in pixels of figure (default 400)
+        :param int width: width in pixels of figure (default 400)
+        """
 
         # Store height and width
         self.width = width
         self.height = height
         # Make temporary folder for all our files
-        self.tmp_dir = Path("./tmp/")
-        try:
-            os.mkdir(self.tmp_dir)
-        except FileExistsError:
-            pass
-        # Create new environment for tmp WGS84 location
-        rcfile, self._vector_env = gs.create_environment(
-            self.tmp_dir, "temp_folium_WGS84", "PERMANENT"
-        )
-        # Location and mapset and region
-        gs.create_location(
-            self.tmp_dir, "temp_folium_WGS84", epsg="4326", overwrite=True
+        self._tmp_dir = tempfile.TemporaryDirectory()
+
+        # Remember original environment; all environments used
+        # in this class are derived from this one
+        self._src_env = os.environ.copy()
+
+        # Set up temporary locations  in WGS84 and Pseudo-Mercator
+        # We need two because folium uses WGS84 for vectors and coordinates
+        # and Pseudo-Mercator for raster overlays
+        self.rcfile_psmerc, self._psmerc_env = setup_location(
+            "psmerc", self._tmp_dir.name, "3857", self._src_env
         )
-        self._extent = self._convert_extent(
-            env=os.environ
-        )  # Get the extent of the original area in WGS84
-        # Set region to match original region extent
-        gs.run_command(
-            "g.region",
-            n=self._extent["north"],
-            s=self._extent["south"],
-            e=self._extent["east"],
-            w=self._extent["west"],
-            env=self._vector_env,
+        self.rcfile_wgs84, self._wgs84_env = setup_location(
+            "wgs84", self._tmp_dir.name, "4326", self._src_env
         )
-        # Get Center of tmp GRASS region
-        center = gs.parse_command("g.region", flags="cg", env=self._vector_env)
+
+        # Get Center of temporary GRASS regions
+        center = gs.parse_command("g.region", flags="cg", env=self._wgs84_env)
         center = (float(center["center_northing"]), float(center["center_easting"]))
+
         # Create Folium Map
         self.map = folium.Map(
             width=self.width,
@@ -67,87 +76,130 @@ class InteractiveMap:
             location=center,
             tiles="cartodbpositron",
         )
-        # Create LayerControl default
+        # Set LayerControl default
         self.layer_control = False
 
-    def _convert_coordinates(self, x, y, proj_in):
-        """This function reprojects coordinates to WGS84, the required
-        projection for vectors in folium.
-
-        Arguments:
-            x -- x coordinate (string)
-            y -- y coordinate (string)
-            proj_in -- proj4 string of location (for example, the output
-            of g.region run with the `g` flag."""
-
-        # Reformat input
-        coordinates = f"{x}, {y}"
-        # Reproject coordinates
-        coords_folium = gs.read_command(
-            "m.proj",
-            coordinates=coordinates,
-            proj_in=proj_in,
-            separator="comma",
-            flags="do",
+        # Cleanup rcfiles with finalizer
+        def remove_if_exists(path):
+            if sys.version_info < (3, 8):
+                try:
+                    os.remove(path)
+                except FileNotFoundError:
+                    pass
+            else:
+                path.unlink(missing_ok=True)
+
+        def clean_up(paths):
+            for path in paths:
+                remove_if_exists(path)
+
+        self._finalizer = weakref.finalize(
+            self, clean_up, [Path(self.rcfile_psmerc), Path(self.rcfile_wgs84)]
         )
-        # Reformat from string to array
-        coords_folium = coords_folium.strip()  # Remove '\n' at end of string
-        coords_folium = coords_folium.split(",")  # Split on comma
-        coords_folium = [float(value) for value in coords_folium]  # Convert to floats
-        return coords_folium[1], coords_folium[0]  # Return Lat and Lon
-
-    def _convert_extent(self, env=None):
-        """This function returns the bounding box of the current region
-        in WGS84, the required projection for folium"""
-
-        # Get proj of current GRASS region
-        proj = gs.read_command("g.proj", flags="jf", env=env)
-        # Get extent
-        extent = gs.parse_command("g.region", flags="g", env=env)
-        # Convert extent to EPSG:3857, required projection for Folium
-        north, east = self._convert_coordinates(extent["e"], extent["n"], proj)
-        south, west = self._convert_coordinates(extent["w"], extent["s"], proj)
-        extent = {"north": north, "south": south, "east": east, "west": west}
-        return extent
-
-    def _folium_bounding_box(self, extent):
-        """Reformats extent into bounding box to pass to folium"""
-        return [[extent["north"], extent["west"]], [extent["south"], extent["east"]]]
 
     def add_vector(self, name):
         """Imports vector into temporary WGS84 location,
         re-formats to a GeoJSON and adds to folium map.
 
-        Arguments:
-            name -- a positional-only parameter; name of vector to be added
-            to map as a string"""
+        :param str name: name of vector to be added to map;
+                         positional-only parameter
+        """
 
         # Find full name of vector
         file_info = gs.find_file(name, element="vector")
         full_name = file_info["fullname"]
         name = file_info["name"]
         # Reproject vector into WGS84 Location
-        env_info = gs.gisenv(env=os.environ)
+        env_info = gs.gisenv(env=self._src_env)
         gs.run_command(
             "v.proj",
             input=full_name,
             location=env_info["LOCATION_NAME"],
             dbase=env_info["GISDBASE"],
-            env=self._vector_env,
+            env=self._wgs84_env,
         )
         # Convert to GeoJSON
-        json_file = self.tmp_dir / f"tmp_{name}.json"
+        json_file = Path(self._tmp_dir.name) / f"tmp_{name}.json"
         gs.run_command(
             "v.out.ogr",
             input=name,
             output=json_file,
             format="GeoJSON",
-            env=self._vector_env,
+            env=self._wgs84_env,
         )
         # Import GeoJSON to folium and add to map
         folium.GeoJson(str(json_file), name=name).add_to(self.map)
 
+    def add_raster(self, name, opacity=0.8):
+        """Imports raster into temporary WGS84 location,
+        exports as png and overlays on folium map
+
+        :param str name: name of raster to add to display; positional-only parameter
+        :param float opacity: raster opacity, number between
+                              0 (transparent) and 1 (opaque)
+        """
+
+        # Find full name of raster
+        file_info = gs.find_file(name, element="cell")
+        full_name = file_info["fullname"]
+        name = file_info["name"]
+
+        # Reproject raster into WGS84/epsg3857 location
+        env_info = gs.gisenv(env=self._src_env)
+        resolution = estimate_resolution(
+            raster=full_name,
+            dbase=env_info["GISDBASE"],
+            location=env_info["LOCATION_NAME"],
+            env=self._psmerc_env,
+        )
+        tgt_name = full_name.replace("@", "_")
+        gs.run_command(
+            "r.proj",
+            input=full_name,
+            output=tgt_name,
+            location=env_info["LOCATION_NAME"],
+            dbase=env_info["GISDBASE"],
+            resolution=resolution,
+            env=self._psmerc_env,
+        )
+        # Write raster to png file with GrassRenderer
+        region_info = gs.region(env=self._src_env)
+        png_width = region_info["cols"]
+        png_height = region_info["rows"]
+        filename = os.path.join(self._tmp_dir.name, f"{tgt_name}.png")
+        m = GrassRenderer(
+            width=png_width,
+            height=png_height,
+            env=self._psmerc_env,
+            filename=filename,
+        )
+        m.run("d.rast", map=tgt_name)
+
+        # Reproject bounds of raster for overlaying png
+        # Bounds need to be in WGS84
+        old_bounds = get_region(self._src_env)
+        from_proj = get_location_proj_string(env=self._src_env)
+        to_proj = get_location_proj_string(env=self._wgs84_env)
+        bounds = reproject_region(old_bounds, from_proj, to_proj)
+        new_bounds = [
+            [bounds["north"], bounds["west"]],
+            [bounds["south"], bounds["east"]],
+        ]
+
+        # Overlay image on folium map
+        img = folium.raster_layers.ImageOverlay(
+            image=filename,
+            name=name,
+            bounds=new_bounds,
+            opacity=opacity,
+            interactive=True,
+            cross_origin=False,
+        )
+        # Add image to map
+        img.add_to(self.map)
+
     def add_layer_control(self, **kwargs):
+        """Add layer control to display"""
         self.layer_control = True
         self.layer_control_object = folium.LayerControl(**kwargs)
 
@@ -166,3 +218,10 @@ class InteractiveMap:
         fig.add_child(self.map)
 
         return fig
+
+    def save(self, filename):
+        """Save map as an html map.
+
+        :param str filename: name of html file
+        """
+        self.map.save(filename)

python/grass/jupyter/utils.py

@@ -0,0 +1,127 @@
+#
+# AUTHOR(S): Caitlin Haedrich <caitlin DOT haedrich AT gmail>
+#
+# PURPOSE:   This module contains utility functions for InteractiveMap.
+#
+# COPYRIGHT: (C) 2021 Caitlin Haedrich, and by the GRASS Development Team
+#
+#            This program is free software under the GNU General Public
+#            License (>=v2). Read the file COPYING that comes with GRASS
+#            for details.
+
+import os
+import grass.script as gs
+
+
+def get_region(env=None):
+    """Returns current computational region as dictionary.
+    Adds long key names.
+    """
+    region = gs.region(env=env)
+    region["east"] = region["e"]
+    region["west"] = region["w"]
+    region["north"] = region["n"]
+    region["south"] = region["s"]
+    return region
+
+
+def get_location_proj_string(env=None):
+    """Returns projection of environment in PROJ.4 format"""
+    out = gs.read_command("g.proj", flags="jf", env=env)
+    return out.strip()
+
+
+def reproject_region(region, from_proj, to_proj):
+    """Reproject boundary of region from one projection to another.
+
+    :param dict region: region to reproject as a dictionary with long key names
+                    output of get_region
+    :param str from_proj: PROJ.4 string of region; output of get_location_proj_string
+    :param str in_proj: PROJ.4 string of target location;
+                    output of get_location_proj_string
+
+    :return dict region: reprojected region as a dictionary with long key names
+    """
+    region = region.copy()
+    proj_input = "{east} {north}\n{west} {south}".format(**region)
+    proc = gs.start_command(
+        "m.proj",
+        input="-",
+        separator=" , ",
+        proj_in=from_proj,
+        proj_out=to_proj,
+        flags="d",
+        stdin=gs.PIPE,
+        stdout=gs.PIPE,
+        stderr=gs.PIPE,
+    )
+    proc.stdin.write(gs.encode(proj_input))
+    proc.stdin.close()
+    proc.stdin = None
+    proj_output, stderr = proc.communicate()
+    if proc.returncode:
+        raise RuntimeError("reprojecting region: m.proj error: " + stderr)
+    enws = gs.decode(proj_output).split(os.linesep)
+    elon, nlat, unused = enws[0].split(" ")
+    wlon, slat, unused = enws[1].split(" ")
+    region["east"] = elon
+    region["north"] = nlat
+    region["west"] = wlon
+    region["south"] = slat
+    return region
+
+
+def estimate_resolution(raster, dbase, location, env):
+    """Estimates resolution of reprojected raster.
+
+    :param str raster: name of raster
+    :param str dbase: path to source database
+    :param str location: name of source location
+    :param dict env: target environment
+
+    :return float estimate: estimated resolution of raster in destination
+                            environment
+    """
+    output = gs.read_command(
+        "r.proj", flags="g", input=raster, dbase=dbase, location=location, env=env
+    ).strip()
+    params = gs.parse_key_val(output, vsep=" ")
+    output = gs.read_command("g.region", flags="ug", env=env, **params)
+    output = gs.parse_key_val(output, val_type=float)
+    cell_ns = (output["n"] - output["s"]) / output["rows"]
+    cell_ew = (output["e"] - output["w"]) / output["cols"]
+    estimate = (cell_ew + cell_ns) / 2.0
+    return estimate
+
+
+def setup_location(name, path, epsg, src_env):
+    """Setup temporary location with different projection but
+    same computational region as source location
+
+    :param str name: name of new location
+    :param path path: path to new location's database
+    :param str epsg: EPSG code
+    :param dict src_env: source environment
+
+    :return str rcfile: name of new locations rcfile
+    :return dict new_env: new environment
+    """
+    # Create new environment
+    rcfile, new_env = gs.create_environment(path, name, "PERMANENT")
+    # Location and mapset
+    gs.create_location(path, name, epsg=epsg, overwrite=True)
+    # Reproject region
+    region = get_region(env=src_env)
+    from_proj = get_location_proj_string(src_env)
+    to_proj = get_location_proj_string(env=new_env)
+    new_region = reproject_region(region, from_proj, to_proj)
+    # Set region to match original region extent
+    gs.run_command(
+        "g.region",
+        n=new_region["north"],
+        s=new_region["south"],
+        e=new_region["east"],
+        w=new_region["west"],
+        env=new_env,
+    )
+    return rcfile, new_env