Testing examples from hydroclimatology#

#!conda env update -f '../../environment.yml' --prune

# this block is for developing the module, comment out when using the module, and uncomment import easysnowdata
%load_ext autoreload
%autoreload 2
%aimport easysnowdata

#import easysnowdata
import geopandas as gpd
import xarray as xr
#from easysnowdata.utils import convert_bbox_to_geodataframe
import json
import ee
import matplotlib.pyplot as plt
import numpy as np
import contextily as ctx
import rioxarray as rxr
import pandas as pd

ee.Authenticate()  # need to figure out https://developers.google.com/earth-engine/guides/auth
ee.Initialize(project='egagli-data-access',opt_url="https://earthengine-highvolume.googleapis.com")

*** Earth Engine *** Share your feedback by taking our Annual Developer Satisfaction Survey: https://google.qualtrics.com/jfe/form/SV_7TDKVSyKvBdmMqW?ref=4i2o6

bbox_gdf = gpd.read_file(
    "https://github.com/egagli/easysnowdata/raw/main/docs/examples/mt_rainier.geojson"
)

Koppen Geiger Climate Classification #

koppen_geiger_da = easysnowdata.hydroclimatology.get_koppen_geiger_classes(resolution="0.1 degree")
#koppen_geiger_da

No spatial subsetting because bbox_input was not provided.

f,ax = koppen_geiger_da.attrs['example_plot'](koppen_geiger_da)

../_images/b1c074260da786cd138cccecfb6113547d957fe272b5a1232b0a0764772eeb61.png

koppen_geiger_da = easysnowdata.hydroclimatology.get_koppen_geiger_classes(bbox_input=bbox_gdf,resolution="1 km")
#koppen_geiger_da

f,ax = koppen_geiger_da.attrs['example_plot'](koppen_geiger_da)

../_images/79a72ae47ef882ad5250b43c9aafee77b39af138d9f620a87bf86973f8f1b669.png

HUC02 geometries #

This might take a moment…

huc02_gdf = easysnowdata.hydroclimatology.get_huc_geometries(huc_level="02")

No spatial subsetting because bbox_input was not provided.

# huc02_gdf.explore()

Get HUC geometries given a bounding box#

huc02_gdf = easysnowdata.hydroclimatology.get_huc_geometries(
    bbox_input=bbox_gdf, huc_level="02"
)
huc04_gdf = easysnowdata.hydroclimatology.get_huc_geometries(
    bbox_input=bbox_gdf, huc_level="04"
)
huc06_gdf = easysnowdata.hydroclimatology.get_huc_geometries(
    bbox_input=bbox_gdf, huc_level="06"
)
huc08_gdf = easysnowdata.hydroclimatology.get_huc_geometries(
    bbox_input=bbox_gdf, huc_level="08"
)
huc10_gdf = easysnowdata.hydroclimatology.get_huc_geometries(
    bbox_input=bbox_gdf, huc_level="10"
)
huc12_gdf = easysnowdata.hydroclimatology.get_huc_geometries(
    bbox_input=bbox_gdf, huc_level="12"
)

f, axs = plt.subplots(1, 2, figsize=(15, 7))

for ax in axs:
    huc02_gdf.plot(ax=ax, color="none", edgecolor="red", linewidth=2, label="HUC02")
    huc04_gdf.plot(ax=ax, color="none", edgecolor="orange", linewidth=2, label="HUC04")
    huc06_gdf.plot(ax=ax, color="none", edgecolor="yellow", linewidth=2, label="HUC06")
    huc08_gdf.plot(ax=ax, color="none", edgecolor="green", linewidth=2, label="HUC08")
    huc10_gdf.plot(ax=ax, color="none", edgecolor="blue", linewidth=2, label="HUC10")
    huc12_gdf.plot(ax=ax, color="none", edgecolor="purple", linewidth=2, label="HUC12")
    bbox_gdf.plot(ax=ax, color="none", edgecolor="black", linewidth=3, label="bbox")

ctx.add_basemap(axs[0], crs=huc02_gdf.crs)

axs[1].set_xlim(-122.5, -121.1)
axs[1].set_ylim(46.4, 47.3)
ctx.add_basemap(axs[1], crs=huc02_gdf.crs)


# Create a legend for the HUC levels
from matplotlib.lines import Line2D

legend_elements = [
    Line2D([0], [0], color="black", lw=2, label="bbox"),
    Line2D([0], [0], color="red", lw=2, label="HUC02"),
    Line2D([0], [0], color="orange", lw=2, label="HUC04"),
    Line2D([0], [0], color="yellow", lw=2, label="HUC06"),
    Line2D([0], [0], color="green", lw=2, label="HUC08"),
    Line2D([0], [0], color="blue", lw=2, label="HUC10"),
    Line2D([0], [0], color="purple", lw=2, label="HUC12"),
]
for ax in axs:
    ax.legend(handles=legend_elements)

f.tight_layout()
f.suptitle("HUC geometries")

Text(0.5, 0.98, 'HUC geometries')

../_images/1859f9d6c8043e8c428246828f3e1b1cedb5f2541cdbacacdb3351980f67d62f.png

HydroBASINS geometries #

hydroATLAS_level5_gdf = easysnowdata.hydroclimatology.get_hydroBASINS(bbox_input=bbox_gdf, level=5)
hydroATLAS_level5_gdf

Loading HydroATLAS level 5 basins...

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/geopandas/io/file.py:497: UserWarning: More than one layer found in 'BasinATLAS_Data_v10.gdb.zip': 'BasinATLAS_v10_lev01' (default), 'BasinATLAS_v10_lev02', 'BasinATLAS_v10_lev03', 'BasinATLAS_v10_lev04', 'BasinATLAS_v10_lev05', 'BasinATLAS_v10_lev06', 'BasinATLAS_v10_lev07', 'BasinATLAS_v10_lev08', 'BasinATLAS_v10_lev09', 'BasinATLAS_v10_lev10', 'BasinATLAS_v10_lev11', 'BasinATLAS_v10_lev12'. Specify layer parameter to avoid this warning.
  crs = pyogrio.read_info(path_or_bytes).get("crs")

	HYBAS_ID	NEXT_DOWN	NEXT_SINK	MAIN_BAS	DIST_SINK	DIST_MAIN	SUB_AREA	UP_AREA	PFAF_ID	ENDO	...	hft_ix_s09	hft_ix_u09	gad_id_smj	gdp_ud_sav	gdp_ud_ssu	gdp_ud_usu	hdi_ix_sav	Shape_Length	Shape_Area	geometry
0	7.050015e+09	0.0	7.050015e+09	7.050015e+09	0.0	0.0	12236.9	651407.2	78210	0	...	63	42	240	54861	2.388950e+10	3.902610e+11	946	7.511620	1.430891	MULTIPOLYGON (((-123.34583 46.09167, -123.3647...
1	7.050015e+09	0.0	7.050015e+09	7.050015e+09	0.0	0.0	52968.8	52968.8	78310	0	...	89	38	240	55345	3.008071e+11	1.575763e+07	947	44.107357	6.367568	MULTIPOLYGON (((-123.96806 46.25833, -123.9951...

2 rows × 297 columns

f, ax = plt.subplots(figsize=(12, 8))
hydroATLAS_level5_gdf.plot(ax=ax, column='HYBAS_ID', color='none', edgecolor='red', linewidth=2)
bbox_gdf.plot(ax=ax, color='none', edgecolor='black', linewidth=3, label='bbox')
ctx.add_basemap(ax, crs=hydroATLAS_level5_gdf.crs)
ax.set_title('HydroATLAS Level 5 Basins - Mt. Rainier Region')
ax.legend()

/tmp/ipykernel_356263/1465880548.py:2: UserWarning: Only specify one of 'column' or 'color'. Using 'color'.
  hydroATLAS_level5_gdf.plot(ax=ax, column='HYBAS_ID', color='none', edgecolor='red', linewidth=2)
/tmp/ipykernel_356263/1465880548.py:6: UserWarning: Legend does not support handles for PatchCollection instances.
See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html#implementing-a-custom-legend-handler
  ax.legend()
/tmp/ipykernel_356263/1465880548.py:6: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  ax.legend()

<matplotlib.legend.Legend at 0x7f1284a0f0a0>

../_images/e20b804b3f05380e23610155fb08557f7806b68810a806cf396c1a4496dc7ba4.png

hydroATLAS_level3_gdf = easysnowdata.hydroclimatology.get_hydroBASINS(bbox_input=bbox_gdf, level=3)
hydroATLAS_level7_gdf = easysnowdata.hydroclimatology.get_hydroBASINS(bbox_input=bbox_gdf, level=7)

f, axs = plt.subplots(1, 2, figsize=(15, 7))

for ax in axs:
    hydroATLAS_level3_gdf.plot(ax=ax, color="none", edgecolor="red", linewidth=2, label="Level 3")
    hydroATLAS_level5_gdf.plot(ax=ax, color="none", edgecolor="orange", linewidth=2, label="Level 5")
    hydroATLAS_level7_gdf.plot(ax=ax, color="none", edgecolor="blue", linewidth=2, label="Level 7")
    bbox_gdf.plot(ax=ax, color="none", edgecolor="black", linewidth=3, label="bbox")

ctx.add_basemap(axs[0], crs=hydroATLAS_level5_gdf.crs)

axs[1].set_xlim(-122.5, -121.1)
axs[1].set_ylim(46.4, 47.3)
ctx.add_basemap(axs[1], crs=hydroATLAS_level5_gdf.crs)

# Create a legend for the HydroATLAS levels
from matplotlib.lines import Line2D

legend_elements = [
    Line2D([0], [0], color="black", lw=3, label="bbox"),
    Line2D([0], [0], color="red", lw=2, label=f"Level 3 ({len(hydroATLAS_level3_gdf)} basins)"),
    Line2D([0], [0], color="orange", lw=2, label=f"Level 5 ({len(hydroATLAS_level5_gdf)} basins)"),
    Line2D([0], [0], color="blue", lw=2, label=f"Level 7 ({len(hydroATLAS_level7_gdf)} basins)"),
]
for ax in axs:
    ax.legend(handles=legend_elements)

f.tight_layout()
f.suptitle("HydroATLAS Basin Hierarchy - Mt. Rainier Region")

Loading HydroATLAS level 3 basins...
Loading HydroATLAS level 7 basins...

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/geopandas/io/file.py:497: UserWarning: More than one layer found in 'BasinATLAS_Data_v10.gdb.zip': 'BasinATLAS_v10_lev01' (default), 'BasinATLAS_v10_lev02', 'BasinATLAS_v10_lev03', 'BasinATLAS_v10_lev04', 'BasinATLAS_v10_lev05', 'BasinATLAS_v10_lev06', 'BasinATLAS_v10_lev07', 'BasinATLAS_v10_lev08', 'BasinATLAS_v10_lev09', 'BasinATLAS_v10_lev10', 'BasinATLAS_v10_lev11', 'BasinATLAS_v10_lev12'. Specify layer parameter to avoid this warning.
  crs = pyogrio.read_info(path_or_bytes).get("crs")
/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/geopandas/io/file.py:497: UserWarning: More than one layer found in 'BasinATLAS_Data_v10.gdb.zip': 'BasinATLAS_v10_lev01' (default), 'BasinATLAS_v10_lev02', 'BasinATLAS_v10_lev03', 'BasinATLAS_v10_lev04', 'BasinATLAS_v10_lev05', 'BasinATLAS_v10_lev06', 'BasinATLAS_v10_lev07', 'BasinATLAS_v10_lev08', 'BasinATLAS_v10_lev09', 'BasinATLAS_v10_lev10', 'BasinATLAS_v10_lev11', 'BasinATLAS_v10_lev12'. Specify layer parameter to avoid this warning.
  crs = pyogrio.read_info(path_or_bytes).get("crs")

Text(0.5, 0.98, 'HydroATLAS Basin Hierarchy - Mt. Rainier Region')

../_images/cd7ca678b2d43b782c8ae0d2cb534a08c29f2ccee16a750530def90146795b97.png

GDRC Major River Basins of the World #

grdc_basins_gdf = easysnowdata.hydroclimatology.get_grdc_major_river_basins_of_the_world()

No spatial subsetting because bbox_input was not provided.

f, ax = plt.subplots(figsize=(12, 7))
grdc_basins_gdf.to_crs(epsg=3857).plot(ax=ax, column='NAME',cmap='tab20')
ctx.add_basemap(ax, source=ctx.providers.CartoDB.DarkMatter)
ax.set_title("GRDC Major River Basins of the World")

Text(0.5, 1.0, 'GRDC Major River Basins of the World')

../_images/1b6b8776eee2c0adb7b8532e375240d93550b604c3084b1601c00f567dfab6fb.png

GRDC WMO Basins and Sub-Basins #

wmo_basins_gdf = easysnowdata.hydroclimatology.get_grdc_wmo_basins()
#wmo_basins_gdf

No spatial subsetting because bbox_input was not provided.

f, ax = plt.subplots(figsize=(12, 7))
wmo_basins_gdf.plot(ax=ax, column='WMOBB_NAME',cmap='tab20')
ctx.add_basemap(ax, crs=wmo_basins_gdf.crs,source=ctx.providers.CartoDB.DarkMatter)
ax.set_title("GRDC WMO Basins and Sub-basins")

Text(0.5, 1.0, 'GRDC WMO Basins and Sub-basins')

../_images/13d9dce1d04cbf94aacdcc69e1b91ea77170657bec26bb012ae17ae7a9b88ee0.png

UCLA snow reanalysis #

snow_reanalysis_da = easysnowdata.hydroclimatology.get_ucla_snow_reanalysis(
    bbox_gdf,
    start_date="2020-10-01",
    end_date="2021-09-30",
    variable="SWE_Post",
    stats="mean",
)

f = snow_reanalysis_da.isel(time=slice(0, 365, 30)).plot.imshow(
    col="time",
    col_wrap=5,
    cmap="Blues",
    vmin=0,
    vmax=3,
)

f.fig.suptitle("UCLA SWE reanalysis")

Text(0.5, 0.98, 'UCLA SWE reanalysis')

../_images/c2e2f7cf459aa14644dcf6017a7b0e595967f59d01405cb925cfc25994ce7474.png