A Map Visualization Example using PyViz Tools in CyberGIS-Jupyter for water¶
1. Objective¶
This notebook is an example of a interactive map visualization which is the high level visualization using PyViz tools as post-processing of environmental modeling. For this example, we used following PyViz tools:
- geopandas (https://geopandas.org/), cartopy (https://scitools.org.uk/cartopy/), geoviews (https://geoviews.org/), and holoviews (https://holoviews.org/)
2. Watershed: Coweeta subwatershed 18, NC, USA¶
The Coweeta Long Term Ecological Research (LTER) Program is based on a collaborative agreement between the University of Georgia Research Foundation in Athens, Georgia, and the USDA Forest Service Coweeta Hydrologic Laboratory in Macon County, North Carolina. This study area is Coweeta subbasin 18 where is the one of the subbasin in Coweeta watershed. (Coweeta: 16 km2, subbasin18: 0.124km2)
3. Create Model Output (RHESSys Patch) as NetCDF format¶
In [1]:
# import Python library
import pandas as pd
import numpy as np
import xarray as xr
from netCDF4 import Dataset, date2num, num2date
from datetime import datetime, timedelta
import geopandas as gpd
import os
In [2]:
# check the current directory
current_path = os.getcwd()
current_path
Out[2]:
In [3]:
# show files in current directory
os.listdir(current_path)
Out[3]:
In [4]:
!unzip -o rhessys_output_patch.daily.zip
In [5]:
# create pandas dataframe from RHESSys patch output
model_output = pd.read_csv(current_path + "/rhessys_output" +'_patch.daily', delimiter=" ")
model_output.head()
Out[5]:
In [6]:
# set shapefile
shapefile = current_path + "/cw18_proj.shp"
# create geo dataframe using geopandas
shapes_df = gpd.read_file(shapefile, driver='ESRI Shapefile')
In [7]:
# set the value of dimension and coordinate for NetCDF
dates = pd.date_range('1/1/2005', '12/30/2008')
shape = (len(dates), len(shapes_df))
dims = ('time', 'patch', )
patchs = shapes_df['gridcode'].tolist()
coords = {'time': dates, 'patch': patchs}
# Create the data structure of RHESSys mode output
cw18 = xr.Dataset(coords=coords)
cw18
Out[7]:
In [8]:
# Check All Model output variables except for date variables
all_var_col = model_output.columns.values
var_col = all_var_col[4:]
var_col
Out[8]:
In [9]:
# set RHESSys output varialbe :Saturation deficit with depth (unit: mm)
variable = 'sat_def_z'
In [10]:
# set RHESSys output varialbe :Saturation deficit with depth (unit: mm)
variable = 'sat_def_z'
for varname in [variable]:
cw18[varname] = xr.DataArray(data=np.full(shape, np.nan), coords=coords, dims=dims, name=varname)
# show the Dimesions and Coordinates in xarray Dataset
cw18
Out[10]:
In [11]:
# check the Dimetions and Coordinates in xarrya Dataset
cw18
Out[11]:
In [12]:
# extract the "sat_def_z" data from model patch output during certain period (2008-01-01~2008-1-30)
model_output[["year", "month", "day", variable]][365*3*1162:365*3*1162+30*1162]
Out[12]:
In [13]:
# get values of variable and insert multi-array model output considering patch and time dimension into xarray Dataset
single_model_output = model_output[variable].values
multi_array_model_output = single_model_output.reshape(cw18.dims['time'], cw18.dims['patch'])
cw18[variable].values[:, :] = multi_array_model_output
In [14]:
# check a variable in xarrya Dataset
cw18
Out[14]:
In [15]:
# create new NetCDF file for "sat_def_z" data of model patch output
cw18.to_netcdf('RHESSys_output_cw18.nc')
cw18
Out[15]:
4. Create Interactive Map Visualization using NetCDF¶
In [16]:
%pylab inline
import cartopy
import geoviews as gv
import holoviews as hv
hv.notebook_extension('bokeh')
4.1 Display Coweeta subbasin 18 with the background map¶
The background map has a limitation of resolution. So if you enlarge the map, the background map will disappear. But it looks cool!!
In [17]:
shapefile = current_path + "/cw18_proj.shp"
shapes = cartopy.io.shapereader.Reader(shapefile)
(gv.tile_sources.StamenTerrainRetina
* gv.Shape.from_shapefile(shapefile, crs=cartopy.crs.PlateCarree()).opts(
style=dict(fill_color='honeydew', line_color='royalblue', alpha=0.8))
* gv.Points([(-83.435, 35.09),(-83.443, 35.050), (-83.450, 35.02)]).opts(style=dict(size=0))
).opts(width=900, height=500)
Out[17]:
4.2 Create Interactive map using a "sat_def_z" variable at certain date¶
In [18]:
def daily_var_map(var, timestep):
# Convert the data from an xarray dataset to a pandas dataframe, This is necessary for geoviews to be able to plot it
out_df = cw18[var].sel(time=timestep).to_dataframe()
# Make sure we have some metadata to join with the shapefile
out_df['gridcode'] = shapes_df['gridcode'].values
# Create the shape plot - some keys:
# - shapes.records() provides the geometry from the shapefiles
# - out_df provides the data
# - value=var chooses which column of the dataframe to use as data for coloring
# - index=['gridcode'] will provide a name on mouse hover
# - on='gridcode' is the key to join the `shapes` with `out_df`
poly = gv.Shape.from_records(shapes.records(), out_df, value=var, index=['gridcode'], on='gridcode')
# Add some options to make things a bit nicer
# - width=700, sets the width, as we expect
# - cmap='plasma' sets the colormap to 'plasma'
# - tools=['hover'] provides information on mouse hover over
# - colorbar=True adds a colorbar
# - alpha=0.7 adds a bit of transparency
poly = poly.opts(width=700, height=600, cmap='plasma', tools=['hover'], colorbar=True, alpha=0.7)
return poly
In [19]:
# check the desired date
cw18.time.values[100]
Out[19]:
In [20]:
%%time
daily_var_map(variable, cw18.time.values[100])
Out[20]:
4.3 Create Interactive map using a "sat_def_z" variable during certain periods¶
In [21]:
# check the desired date
cw18_desire = [cw18.time.values[365], cw18.time.values[424], cw18.time.values[485], cw18.time.values[546], cw18.time.values[607]]
cw18_desire
Out[21]:
In [22]:
%%time
mapping = {t: daily_var_map('sat_def_z', t) for t in cw18_desire}
hv.HoloMap(mapping, kdims=['timestep'])
Out[22]:
4.4 Create Interactive map using a "sat_def_z" variable with time series graph¶
In [23]:
def interactive_plot(var, start, stop, step):
# Just as before, use the `var_map` function to build our map
poly = hv.HoloMap({t: daily_var_map(var, t) for t in cw18.time.values[start:stop:step]}, kdims=['timestep'])
# Vlines will give us an indicator of which time slice we are looking at
vlines = hv.HoloMap({t: hv.VLine(t) for t in cw18.time.values[start:stop:step]}, kdims=['timestep'])
# This seems to make the plot appear more often - probably a holo/geoviews bug somewhere
poly
# Calculate min, max, and mean over the domain for each time
vmin = cw18[var].isel(time=slice(start, stop)).min(dim='patch').rolling(time=step).min().values
vmax = cw18[var].isel(time=slice(start, stop)).max(dim='patch').rolling(time=step).max().values
vmean = cw18[var].isel(time=slice(start, stop)).mean('patch').rolling(center=True, time=step).mean().dropna('time')
# Build the complete interactive plot. Layers are as follows
# - gv.tile_sources.EsriTerrain- Add a background with topography
# - * poly - overlay our map onto the background
# - + hv.Area... - add the area plot to the right
# - * hv.Curve... - overlay the mean curve onto the area plot
# - * vlines... - overlay the timestep indicator
return (gv.tile_sources.EsriTerrain
* poly
+ (hv
.Area((cw18.time[start:stop], vmin, vmax), vdims=['vmin', 'vmax'])
.opts(alpha=0.5, color='gold', line_color=None)
.redim.label(x='Date', vmin=var.capitalize()))
* hv.Curve(vmean).opts(color='purple', alpha=0.8)
* vlines.opts(alpha=0.4, color='red'))
In [24]:
%%time
interactive_plot(variable, start=365, stop=700, step=60)
Out[24]:
In [ ]: