Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

grass.jupyter: Use InteractiveMap with custom projections #4204

Draft
wants to merge 8 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 37 additions & 12 deletions python/grass/jupyter/interactivemap.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#
# AUTHOR(S): Caitlin Haedrich <caitlin DOT haedrich AT gmail>
# Anna Petrasova <kratochanna AT gmail>
# Riya <29riyasaxena AT gmail>
#
# PURPOSE: This module contains functions for interactive visualizations
# in Jupyter Notebooks.
Expand All @@ -16,7 +17,8 @@
import base64
import json
from pathlib import Path
from .reprojection_renderer import ReprojectionRenderer
from pyproj import CRS
from .reprojection_renderer import ReprojectionRenderer, DirectRenderer

from .utils import (
get_region_bounds_latlon,
Expand Down Expand Up @@ -213,6 +215,7 @@ def __init__(
use_region=False,
saved_region=None,
map_backend=None,
crs=None,
):
"""Creates a blank folium/ipyleaflet map centered on g.region.

Expand Down Expand Up @@ -305,15 +308,41 @@ def _import_ipyleaflet(error):
# Store Region
self.region = None

if self._ipyleaflet:
basemap = xyzservices.providers.query_name(tiles)
if API_key and basemap.get("accessToken"):
basemap["accessToken"] = API_key
layout = self._ipywidgets.Layout(width=f"{width}px", height=f"{height}px")
self.map = self._ipyleaflet.Map(
basemap=basemap, layout=layout, scroll_wheel_zoom=True
if crs:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need to get the CRS of the current project with g.proj, use the 'srid' key to get EPSG

if isinstance(crs, str):
crs_n = CRS.from_string(crs)
crs_dict = {"name": crs_n.to_epsg(), "custom": False}
else:
crs_dict = crs

self._renderer = DirectRenderer(
crs=crs_dict, use_region=use_region, saved_region=saved_region
)

else:
self._renderer = ReprojectionRenderer(
use_region=use_region, saved_region=saved_region
)

if self._ipyleaflet:
if crs:
layout = self._ipywidgets.Layout(
width=f"{width}px", height=f"{height}px"
)
self.map = self._ipyleaflet.Map(
crs=crs_dict, layout=layout, scroll_wheel_zoom=True
)
else:
basemap = xyzservices.providers.query_name(tiles)
if API_key and basemap.get("accessToken"):
basemap["accessToken"] = API_key
layout = self._ipywidgets.Layout(
width=f"{width}px", height=f"{height}px"
)
self.map = self._ipyleaflet.Map(
basemap=basemap, layout=layout, scroll_wheel_zoom=True
)

else:
self.map = self._folium.Map(
width=self.width,
Expand All @@ -325,10 +354,6 @@ def _import_ipyleaflet(error):
self.layer_control_object = None
self.region_rectangle = None

self._renderer = ReprojectionRenderer(
use_region=use_region, saved_region=saved_region
)

def add_vector(self, name, title=None, **kwargs):
"""Imports vector into temporary WGS84 location, re-formats to a GeoJSON and
adds to map.
Expand Down
47 changes: 47 additions & 0 deletions python/grass/jupyter/region.py
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,53 @@ def _set_bbox(self, env):
self._bbox[1][1] = max(self._bbox[1][1], east)


class RegionManagerForDirectRenderer:

def __init__(self, use_region, saved_region, src_env):
self._use_region = use_region
self._saved_region = saved_region
self._src_env = src_env
# [SW, NE]: inverted to easily expand based on data, see _set_bbox
self._bbox = [[90, 180], [-90, -180]]

if self._use_region:
self._set_bbox(self._src_env)

if self._saved_region:
self._src_env["GRASS_REGION"] = gs.region_env(
region=self._saved_region, env=self._src_env
)
self._set_bbox(self._src_env)

@property
def bbox(self):
return self._bbox

def set_region_from_raster(self, raster):
if not (self._use_region or self._saved_region):
self._src_env["GRASS_REGION"] = gs.region_env(
raster=raster, env=self._src_env
)
self._set_bbox(self._src_env)

def set_bbox_vector(self, vector):
if not (self._saved_region or self._use_region):
env = self._src_env.copy()
env["GRASS_REGION"] = gs.region_env(vector=vector, env=env)
self._set_bbox(env)

def _set_bbox(self, env):
bbox = gs.parse_command("g.region", flags="bg", env=env)
south = float(bbox["ll_s"])
west = float(bbox["ll_w"])
north = float(bbox["ll_n"])
east = float(bbox["ll_e"])
self._bbox[0][0] = min(self._bbox[0][0], south)
self._bbox[0][1] = min(self._bbox[0][1], west)
self._bbox[1][0] = max(self._bbox[1][0], north)
self._bbox[1][1] = max(self._bbox[1][1], east)


class RegionManagerFor2D:
"""Region manager for 2D displays (gets region from display commands)"""

Expand Down
69 changes: 66 additions & 3 deletions python/grass/jupyter/reprojection_renderer.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
import tempfile
import weakref
from pathlib import Path

import grass.script as gs
from .map import Map
from .utils import (
Expand All @@ -26,7 +27,7 @@
reproject_region,
setup_location,
)
from .region import RegionManagerForInteractiveMap
from .region import RegionManagerForInteractiveMap, RegionManagerForDirectRenderer


class ReprojectionRenderer:
Expand Down Expand Up @@ -68,10 +69,10 @@ def cleanup(tmpdir):
# 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
"psmerc", self._tmp_dir.name, self._src_env, epsg="3857"
)
self._rcfile_wgs84, self._wgs84_env = setup_location(
"wgs84", self._tmp_dir.name, "4326", self._src_env
"wgs84", self._tmp_dir.name, self._src_env, epsg="4326"
)

# region handling
Expand Down Expand Up @@ -168,3 +169,65 @@ def render_vector(self, name):
)

return json_file


class DirectRenderer:
def __init__(self, crs, use_region=False, saved_region=None, work_dir=None):

if not work_dir:
self._tmp_dir = tempfile.TemporaryDirectory()
else:
self._tmp_dir = work_dir

self._src_env = os.environ.copy()
self._proj_string = get_location_proj_string(env=self._src_env)
self._crs = crs
self._region_manager = RegionManagerForDirectRenderer(
use_region, saved_region, self._src_env
)

def get_bbox(self):
return self._region_manager.bbox

def render_raster(self, name):
file_info = gs.find_file(name, element="cell", env=self._src_env)
full_name = file_info["fullname"]

self._region_manager.set_region_from_raster(full_name)
raster_info = gs.raster_info(full_name, env=self._src_env)
filename = os.path.join(
self._tmp_dir.name, f"{full_name.replace('@', '_')}.png"
)
img = Map(
width=int(raster_info["cols"]),
height=int(raster_info["rows"]),
env=self._src_env,
filename=filename,
use_region=True,
)
img.run("d.rast", map=full_name)

bounds = get_region(self._src_env)
new_bounds = [
[bounds["north"], bounds["west"]],
[bounds["south"], bounds["east"]],
]

return filename, new_bounds

def render_vector(self, name):
file_info = gs.find_file(name, element="vector")
full_name = file_info["fullname"]

self._region_manager.set_bbox_vector(full_name)

json_file = Path(self._tmp_dir.name) / f"{full_name.replace('@', '_')}.json"
gs.run_command(
"v.out.ogr",
input=full_name,
output=json_file,
format="GeoJSON",
env=self._src_env,
)

return json_file
7 changes: 5 additions & 2 deletions python/grass/jupyter/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -331,7 +331,7 @@ def estimate_resolution(raster, mapset, location, dbase, env):
return (cell_ew + cell_ns) / 2.0


def setup_location(name, path, epsg, src_env):
def setup_location(name, path, src_env, epsg=None, proj4=None):
"""Setup temporary location with different projection but
same computational region as source location

Expand All @@ -346,7 +346,10 @@ def setup_location(name, path, epsg, src_env):
# Create new environment
rcfile, new_env = gs.create_environment(path, name, "PERMANENT")
# Location and mapset
gs.create_location(path, name, epsg=epsg, overwrite=True)
if proj4:
gs.create_location(path, name, proj4=proj4, overwrite=True)
elif epsg:
gs.create_location(path, name, epsg=epsg, overwrite=True)
# Reproject region
set_target_region(src_env, new_env)
return rcfile, new_env
Expand Down