interact_display.py 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228
  1. #
  2. # AUTHOR(S): Caitlin Haedrich <caitlin DOT haedrich AT gmail>
  3. #
  4. # PURPOSE: This module contains functions for interactive display
  5. # in Jupyter Notebooks.
  6. #
  7. # COPYRIGHT: (C) 2021 Caitlin Haedrich, and by the GRASS Development Team
  8. #
  9. # This program is free software under the GNU General Public
  10. # License (>=v2). Read the file COPYING that comes with GRASS
  11. # for details.
  12. import os
  13. import sys
  14. import tempfile
  15. import weakref
  16. from pathlib import Path
  17. import folium
  18. import grass.script as gs
  19. from .display import GrassRenderer
  20. from .utils import (
  21. estimate_resolution,
  22. get_location_proj_string,
  23. get_region,
  24. reproject_region,
  25. setup_location,
  26. )
  27. class InteractiveMap:
  28. """This class creates interative GRASS maps with folium.
  29. Basic Usage:
  30. >>> m = InteractiveMap()
  31. >>> m.add_vector("streams")
  32. >>> m.add_raster("elevation")
  33. >>> m.add_layer_control()
  34. >>> m.show()
  35. """
  36. def __init__(self, width=400, height=400):
  37. """Creates a blank folium map centered on g.region.
  38. :param int height: height in pixels of figure (default 400)
  39. :param int width: width in pixels of figure (default 400)
  40. """
  41. # Store height and width
  42. self.width = width
  43. self.height = height
  44. # Make temporary folder for all our files
  45. self._tmp_dir = tempfile.TemporaryDirectory()
  46. # Remember original environment; all environments used
  47. # in this class are derived from this one
  48. self._src_env = os.environ.copy()
  49. # Set up temporary locations in WGS84 and Pseudo-Mercator
  50. # We need two because folium uses WGS84 for vectors and coordinates
  51. # and Pseudo-Mercator for raster overlays
  52. self.rcfile_psmerc, self._psmerc_env = setup_location(
  53. "psmerc", self._tmp_dir.name, "3857", self._src_env
  54. )
  55. self.rcfile_wgs84, self._wgs84_env = setup_location(
  56. "wgs84", self._tmp_dir.name, "4326", self._src_env
  57. )
  58. # Get Center of temporary GRASS regions
  59. center = gs.parse_command("g.region", flags="cg", env=self._wgs84_env)
  60. center = (float(center["center_northing"]), float(center["center_easting"]))
  61. # Create Folium Map
  62. self.map = folium.Map(
  63. width=self.width,
  64. height=self.height,
  65. location=center,
  66. tiles="cartodbpositron",
  67. )
  68. # Set LayerControl default
  69. self.layer_control = False
  70. # Cleanup rcfiles with finalizer
  71. def remove_if_exists(path):
  72. if sys.version_info < (3, 8):
  73. try:
  74. os.remove(path)
  75. except FileNotFoundError:
  76. pass
  77. else:
  78. path.unlink(missing_ok=True)
  79. def clean_up(paths):
  80. for path in paths:
  81. remove_if_exists(path)
  82. self._finalizer = weakref.finalize(
  83. self, clean_up, [Path(self.rcfile_psmerc), Path(self.rcfile_wgs84)]
  84. )
  85. def add_vector(self, name):
  86. """Imports vector into temporary WGS84 location,
  87. re-formats to a GeoJSON and adds to folium map.
  88. :param str name: name of vector to be added to map;
  89. positional-only parameter
  90. """
  91. # Find full name of vector
  92. file_info = gs.find_file(name, element="vector")
  93. full_name = file_info["fullname"]
  94. name = file_info["name"]
  95. # Reproject vector into WGS84 Location
  96. env_info = gs.gisenv(env=self._src_env)
  97. gs.run_command(
  98. "v.proj",
  99. input=full_name,
  100. location=env_info["LOCATION_NAME"],
  101. dbase=env_info["GISDBASE"],
  102. env=self._wgs84_env,
  103. )
  104. # Convert to GeoJSON
  105. json_file = Path(self._tmp_dir.name) / f"tmp_{name}.json"
  106. gs.run_command(
  107. "v.out.ogr",
  108. input=name,
  109. output=json_file,
  110. format="GeoJSON",
  111. env=self._wgs84_env,
  112. )
  113. # Import GeoJSON to folium and add to map
  114. folium.GeoJson(str(json_file), name=name).add_to(self.map)
  115. def add_raster(self, name, opacity=0.8):
  116. """Imports raster into temporary WGS84 location,
  117. exports as png and overlays on folium map
  118. :param str name: name of raster to add to display; positional-only parameter
  119. :param float opacity: raster opacity, number between
  120. 0 (transparent) and 1 (opaque)
  121. """
  122. # Find full name of raster
  123. file_info = gs.find_file(name, element="cell")
  124. full_name = file_info["fullname"]
  125. name = file_info["name"]
  126. # Reproject raster into WGS84/epsg3857 location
  127. env_info = gs.gisenv(env=self._src_env)
  128. resolution = estimate_resolution(
  129. raster=full_name,
  130. dbase=env_info["GISDBASE"],
  131. location=env_info["LOCATION_NAME"],
  132. env=self._psmerc_env,
  133. )
  134. tgt_name = full_name.replace("@", "_")
  135. gs.run_command(
  136. "r.proj",
  137. input=full_name,
  138. output=tgt_name,
  139. location=env_info["LOCATION_NAME"],
  140. dbase=env_info["GISDBASE"],
  141. resolution=resolution,
  142. env=self._psmerc_env,
  143. )
  144. # Write raster to png file with GrassRenderer
  145. region_info = gs.region(env=self._src_env)
  146. png_width = region_info["cols"]
  147. png_height = region_info["rows"]
  148. filename = os.path.join(self._tmp_dir.name, f"{tgt_name}.png")
  149. m = GrassRenderer(
  150. width=png_width,
  151. height=png_height,
  152. env=self._psmerc_env,
  153. filename=filename,
  154. )
  155. m.run("d.rast", map=tgt_name)
  156. # Reproject bounds of raster for overlaying png
  157. # Bounds need to be in WGS84
  158. old_bounds = get_region(self._src_env)
  159. from_proj = get_location_proj_string(env=self._src_env)
  160. to_proj = get_location_proj_string(env=self._wgs84_env)
  161. bounds = reproject_region(old_bounds, from_proj, to_proj)
  162. new_bounds = [
  163. [bounds["north"], bounds["west"]],
  164. [bounds["south"], bounds["east"]],
  165. ]
  166. # Overlay image on folium map
  167. img = folium.raster_layers.ImageOverlay(
  168. image=filename,
  169. name=name,
  170. bounds=new_bounds,
  171. opacity=opacity,
  172. interactive=True,
  173. cross_origin=False,
  174. )
  175. # Add image to map
  176. img.add_to(self.map)
  177. def add_layer_control(self, **kwargs):
  178. """Add layer control to display"""
  179. self.layer_control = True
  180. self.layer_control_object = folium.LayerControl(**kwargs)
  181. def show(self):
  182. """This function creates a folium map with a GRASS raster
  183. overlayed on a basemap.
  184. If map has layer control enabled, additional layers cannot be
  185. added after calling show()."""
  186. if self.layer_control:
  187. self.map.add_child(self.layer_control_object)
  188. # Create Figure
  189. fig = folium.Figure(width=self.width, height=self.height)
  190. # Add map to figure
  191. fig.add_child(self.map)
  192. return fig
  193. def save(self, filename):
  194. """Save map as an html map.
  195. :param str filename: name of html file
  196. """
  197. self.map.save(filename)