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:

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]:
'/home/jovyan/work/Downloads/37c0e375e3934fe1a61e1bbc2d05f836/37c0e375e3934fe1a61e1bbc2d05f836/data/contents/2_map_visualization'
In [3]:
# show files in current directory
os.listdir(current_path)
Out[3]:
['rhessys_output_patch.daily',
 'rhessys_output_patch.daily.zip',
 'cw18_proj_resmap.xml',
 'cw18_proj.shp.xml',
 'cw18_proj.sbx',
 '.ipynb_checkpoints',
 'cw18_proj_meta.xml',
 'cw18_proj.shx',
 'cw18_proj.prj',
 'cw18_proj.dbf',
 'background_coweeta_sub18_map.png',
 '__MACOSX',
 'RHESSys_output_cw18.nc',
 'cybergis.png',
 'cw18_proj.shp',
 'cw18_proj.sbn',
 'Data_Processing_example_2_Map_Visualization_Example_using_PyViz_tools_in_CyberGIS_Jupyter_for_water.ipynb']
In [4]:
!unzip -o rhessys_output_patch.daily.zip
Archive:  rhessys_output_patch.daily.zip
  inflating: rhessys_output_patch.daily  
  inflating: __MACOSX/._rhessys_output_patch.daily  
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]:
year month day patchID subsurfaceQnet surfaceQnet detention stormdrainYield return rain_thr ... rtz_fc unsat_storage top12cm_storage top12cm_potential_sat rootdepth soildepth top30cm_storage top30cm_potential_sat top60cm_storage top60cm_potential_sat
0 2005 1 1 42534 144.802905 0.0 0.0 0.0 0.0 0.0 ... 286.406196 0.0 17.064077 54.119188 1000.0 18484.988075 67.711251 135.294926 136.018913 270.579706
1 2005 1 1 41879 131.158183 0.0 0.0 0.0 0.0 0.0 ... 261.455892 0.0 17.081344 54.119188 1000.0 18484.988075 67.950619 135.294926 133.540872 270.579706
2 2005 1 1 43190 179.810526 0.0 0.0 0.0 0.0 0.0 ... 347.711989 0.0 17.014028 54.119188 1000.0 18484.988075 66.372640 135.294926 134.469385 270.579706
3 2005 1 1 41225 92.058033 0.0 0.0 0.0 0.0 0.0 ... 185.305376 0.0 17.251842 54.119188 1000.0 18484.988075 68.081485 135.294926 118.207424 270.579706
4 2005 1 1 40570 125.571288 0.0 0.0 0.0 0.0 0.0 ... 251.428443 0.0 17.081318 54.119188 1000.0 18484.988075 67.847537 135.294926 131.547376 270.579706

5 rows × 31 columns

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]:
<xarray.Dataset>
Dimensions:  (patch: 1162, time: 1460)
Coordinates:
  * time     (time) datetime64[ns] 2005-01-01 2005-01-02 ... 2008-12-30
  * patch    (patch) int64 658 659 1311 1312 1313 ... 42536 42537 43190 43191
Data variables:
    *empty*
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]:
array(['subsurfaceQnet', 'surfaceQnet', 'detention', 'stormdrainYield',
       'return', 'rain_thr', 'thr_recharge', 'cap_drain', 'sat_def_z',
       'sat_def', 'rtzStorage', 'ET', 'treeLAI', 'nontreeLAI',
       'SmartIrrigation', 'rtz_totalvol', 'unsat_fc', 'rtz_fc',
       'unsat_storage', 'top12cm_storage', 'top12cm_potential_sat',
       'rootdepth', 'soildepth', 'top30cm_storage',
       'top30cm_potential_sat', 'top60cm_storage',
       'top60cm_potential_sat'], dtype=object)
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]:
<xarray.Dataset>
Dimensions:    (patch: 1162, time: 1460)
Coordinates:
  * time       (time) datetime64[ns] 2005-01-01 2005-01-02 ... 2008-12-30
  * patch      (patch) int64 658 659 1311 1312 1313 ... 42536 42537 43190 43191
Data variables:
    sat_def_z  (time, patch) float64 nan nan nan nan nan ... nan nan nan nan nan
In [11]:
# check the Dimetions and Coordinates in xarrya Dataset
cw18
Out[11]:
<xarray.Dataset>
Dimensions:    (patch: 1162, time: 1460)
Coordinates:
  * time       (time) datetime64[ns] 2005-01-01 2005-01-02 ... 2008-12-30
  * patch      (patch) int64 658 659 1311 1312 1313 ... 42536 42537 43190 43191
Data variables:
    sat_def_z  (time, patch) float64 nan nan nan nan nan ... nan nan nan nan nan
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]:
year month day sat_def_z
1272390 2008 1 1 5783.071944
1272391 2008 1 1 5295.056387
1272392 2008 1 1 6380.639875
1272393 2008 1 1 5331.593301
1272394 2008 1 1 5754.685782
... ... ... ... ...
1307245 2008 1 30 171.604766
1307246 2008 1 30 0.000000
1307247 2008 1 30 0.000000
1307248 2008 1 30 0.000000
1307249 2008 1 30 0.000000

34860 rows × 4 columns

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]:
<xarray.Dataset>
Dimensions:    (patch: 1162, time: 1460)
Coordinates:
  * time       (time) datetime64[ns] 2005-01-01 2005-01-02 ... 2008-12-30
  * patch      (patch) int64 658 659 1311 1312 1313 ... 42536 42537 43190 43191
Data variables:
    sat_def_z  (time, patch) float64 642.5 583.6 791.8 410.9 ... 0.0 0.0 0.0 0.0
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]:
<xarray.Dataset>
Dimensions:    (patch: 1162, time: 1460)
Coordinates:
  * time       (time) datetime64[ns] 2005-01-01 2005-01-02 ... 2008-12-30
  * patch      (patch) int64 658 659 1311 1312 1313 ... 42536 42537 43190 43191
Data variables:
    sat_def_z  (time, patch) float64 642.5 583.6 791.8 410.9 ... 0.0 0.0 0.0 0.0

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')
Populating the interactive namespace from numpy and matplotlib
/data/cigi/cjw-easybuild/conda/pyrhessys-2021-06/lib/python3.7/site-packages/IPython/core/magics/pylab.py:160: UserWarning: pylab import has clobbered these variables: ['num2date', 'date2num', 'shape', 'datetime']
`%matplotlib` prevents importing * from pylab and numpy
  "\n`%matplotlib` prevents importing * from pylab and numpy"