#!/usr/bin/env python
# coding: utf-8

# # State-level Overview of Sewersheds
# 
# In this Jupyter notebook for Python, we visualize sewersheds at the state level.
# 
# ## Data requirements
# 
# - Approximate positions of sewersheds represented as points
#   * Each point needs to have an attribute with size of the population the given sewershed serves.
# - State boundary

# ## Software setup 
# 
# We will use a couple of standard Python packages and GRASS GIS.
# 
# For now, the setup here assumes Linux. Instructions for Windows are available at [GRASS GIS Jupyter notebooks wiki page](https://grasswiki.osgeo.org/wiki/GRASS_GIS_Jupyter_notebooks#Running_a_Jupyter_notebook_locally).

# In[ ]:


# Import Python standard library and IPython packages we need.
import os
import subprocess
import sys
from pathlib import Path

# Ask GRASS GIS where its Python packages are.
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj


# ## Get the data ready
# 
# This notebooks needs sewershed locations as vector points with at attributes for population size and name. The file should be in directory `data/sewers_state` and should be named `sewer_points` with file extension appropriate for the format, e.g. `sewer_points.shp`. Either rename the files or modify the code below if needed.

# In[ ]:


data_directory = "data/sewers_state"
point_file = "sewer_points"  # Filename without file extension
state_boundary_file = "state_boundary"


# In[ ]:


grass_project = "data/state_overview"


# To compute the data, we will use a GRASS project (aka location).

# In[ ]:


get_ipython().system('grass -e -c $data_directory $grass_project')


# In[ ]:


session = gj.init(grass_project)


# From now on, we will use the following variables to refer to the data in the GRASS project.

# In[ ]:


point_vector = "points"
state_boundary_vector = "state_boundary_vector"


# In[ ]:


gs.run_command("v.import", input=data_directory, layer=point_file, output=point_vector)


# In[ ]:


gs.run_command(
    "v.import",
    input=data_directory,
    layer=state_boundary_file,
    output=state_boundary_vector,
)


# ## Render

# In[ ]:


# Extract the large sewersheds.
large_sewersheds = "large_sewersheds"

gs.run_command(
    "v.extract",
    input=point_vector,
    output=large_sewersheds,
    where="pop_2020 > 100000",
)


# In[ ]:


gs.run_command("g.region", vector="nc_state", res=1000, grow=50)
state_map = gj.Map(width=1024, use_region=True)
state_map.d_background(color="white")
state_map.d_vect(flags="s", map=state_boundary_vector, fill_color="none")
state_map.d_vect_thematic(
    map=point_vector,
    column="pop_2020",
    legend_title="Population per sewershed",
    colors="#ADA8A7,#6A00A8,#B12A90,#E16462,#FCA636,#F0F921",
    breaks=[100000, 200000, 300000, 400000],
    boundary_color="none",
    icon="basic/circle",
    size=10,
)
# We rendered the points colored by population. Now we take the large
# sewersheds and render their names as labels in the map (this happens in two
# steps (first, labels are extracted and laid out, then rendered in
# the second step).
gs.run_command(
    "v.label.sa",
    map=large_sewersheds,
    column="site_name",
    labels="points_site_name",
    font="DejaVuSans",
    size=30000,
    color="black",
    isize=25000,
)
state_map.d_labels(labels="points_site_name")
state_map.d_legend_vect(flags="b", at=(3, 34), columns=2)
state_map.d_barscale(flags="n", at=(2, 7), length=200, units="kilometers")
state_map.save(Path(data_directory) / "nc_sewersheds.png")
state_map.show()