Testing examples from hydroclimatology¶

In [1]:

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

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

In [2]:

#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 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

In [3]:

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

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

In [4]:

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

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

Koppen Geiger Climate Classification¶

In [5]:

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

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.

In [6]:

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

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

In [7]:

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

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

In [8]:

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

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

HUC02 geometries¶

This might take a moment...

In [9]:

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

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

No spatial subsetting because bbox_input was not provided.

In [10]:

# huc02_gdf.explore()

# huc02_gdf.explore()

Get HUC geometries given a bounding box¶

In [11]:

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"
)

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" )

In [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")

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")

ERA5 atmospheric reanalysis¶

In [13]:

era5_global_ds = easysnowdata.hydroclimatology.get_era5()

era5_global_ds = easysnowdata.hydroclimatology.get_era5()

No spatial subsetting because bbox_input was not provided.

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/xarray/core/indexing.py:1620: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]

In [14]:

f, ax = plt.subplots()
era5_global_ds["2m_temperature"].sel(time="2020-05-26").mean(dim="time").plot.imshow(ax=ax)
f.suptitle("ERA5 2m Temperature on 2020-05-26")

f, ax = plt.subplots() era5_global_ds["2m_temperature"].sel(time="2020-05-26").mean(dim="time").plot.imshow(ax=ax) f.suptitle("ERA5 2m Temperature on 2020-05-26")

Out[14]:

<matplotlib.image.AxesImage at 0x7fee612862f0>

In [15]:

era5_clipped_ds = easysnowdata.hydroclimatology.get_era5(bbox_gdf)

era5_clipped_ds = easysnowdata.hydroclimatology.get_era5(bbox_gdf)

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/xarray/core/indexing.py:1620: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]

In [16]:

f, ax = plt.subplots()
era5_clipped_ds["2m_temperature"].isel(time=-1).plot.imshow(ax=ax)
bbox_gdf.plot(ax=ax, color="none", edgecolor="red", linewidth=3)

f, ax = plt.subplots() era5_clipped_ds["2m_temperature"].isel(time=-1).plot.imshow(ax=ax) bbox_gdf.plot(ax=ax, color="none", edgecolor="red", linewidth=3)

Out[16]:

<Axes: title={'center': 'time = 2021-12-31T23:00:00, spatial_ref = 0'}, xlabel='longitude [degrees_east]', ylabel='latitude [degrees_north]'>

UCLA snow reanalysis¶

In [17]:

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",
)

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", )

Granules found: 2
Opening 2 granules, approx size: 0.33 GB

QUEUEING TASKS | :   0%|          | 0/2 [00:00<?, ?it/s]

PROCESSING TASKS | :   0%|          | 0/2 [00:00<?, ?it/s]

COLLECTING RESULTS | :   0%|          | 0/2 [00:00<?, ?it/s]

In [18]:

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

f.suptitle("UCLA SWE reanalysis")

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

Out[18]:

<xarray.plot.facetgrid.FacetGrid at 0x7fee40b839a0>

In [ ]: