Basin stats examples
Generating basin stats from Earth Engine imagery is challenged by the many different ways in which image assets are stored and organized on the platform. As such, we have tried to strike a balance with rabpro between ease of use (abstracting away unncessary details) and flexibility. Before we begin pulling basin stats, lets create a basin polygon asset:
[2]:
import pandas as pd
import geopandas as gpd
import rabpro
from rabpro.basin_stats import Dataset
[3]:
# # Authenticate earthengine if necessary
import ee
# ee.Authenticate()
ee.Initialize()
[4]:
%%capture
coords = (44.9331, -69.4996)
rpo = rabpro.profiler(coords, name='basic_test')
rpo.delineate_basin()
gdf = rpo.watershed
The next several examples pull from imagery assets that are within the standard GEE data catalog (https://developers.google.com/earth-engine/datasets/). This simplifies their queries as we (typically) do not need to specify projection or resolution information. These are read directly from the catalog:
Categorical asset - a single time window
[4]:
urls, tasks = rabpro.basin_stats.compute(
[Dataset("MODIS/006/MCD12Q1", "LC_Type1", stats=["freqhist"], start="2010-01-01", end="2011-01-01")],
validate_dataset_list=False,
gee_feature_path="users/jstacompute/rpo_basic"
)
data = rabpro.basin_stats.fetch_gee(urls, ["lulc"], ["system:index"])
res = gpd.GeoDataFrame(pd.concat([data, gdf], axis=1))
res.head()
[4]:
| lulc_1 | lulc_17 | lulc_4 | lulc_5 | lulc_8 | da_km2 | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | 1.87451 | 16.337255 | 621.25098 | 1143.831373 | 262.792157 | 440.266532 | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
Numeric asset - multiple time windows
[5]:
urls, tasks = rabpro.basin_stats.compute(
[Dataset("MODIS/006/MOD17A3HGF", "Npp")],
validate_dataset_list=False,
gee_feature_path="users/jstacompute/rpo_basic"
)
data = rabpro.basin_stats.fetch_gee(urls, ["npp"], ["da_km2"])
[6]:
data["year"] = [x for x in range(2000, 2000 + data.shape[0])]
res = gpd.GeoDataFrame(pd.concat([data, gdf], axis=1))
res.geometry[res.geometry.isna()] = res.geometry[0]
res.head()
[6]:
| npp_system:index | npp_mean | year | da_km2 | geometry | |
|---|---|---|---|---|---|
| 0 | 2000_02_18_00000000000000000000 | 7436.137468 | 2000 | 440.266532 | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
| 1 | 2001_01_01_00000000000000000000 | 7253.867450 | 2001 | NaN | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
| 2 | 2002_01_01_00000000000000000000 | 6741.378228 | 2002 | NaN | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
| 3 | 2003_01_01_00000000000000000000 | 6565.773486 | 2003 | NaN | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
| 4 | 2004_01_01_00000000000000000000 | 7645.184775 | 2004 | NaN | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
Numeric asset - time averaging
[7]:
urls, tasks = rabpro.basin_stats.compute(
[Dataset("MODIS/061/MOD17A3HGF", "Npp", time_stats=["median"])],
gee_feature_path="users/jstacompute/rpo_basic"
)
data = rabpro.basin_stats.fetch_gee(urls, ["npp"], ["system:index"])
res = gpd.GeoDataFrame(pd.concat([data, gdf], axis=1))
res.head()
[7]:
| npp_mean | da_km2 | geometry | |
|---|---|---|---|
| 0 | 7238.642484 | 440.266532 | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |
The final example here demonstates a query of a “custom” imagery asset which is not present in the GEE data catalog. As a result, we must specify additional information to the compute function:
[7]:
urls, tasks = rabpro.basin_stats.compute(
[Dataset("projects/soilgrids-isric/soc_mean", "soc_0-5cm_mean", gee_type="image")],
validate_dataset_list=False,
gee_feature_path="users/jstacompute/rpo_basic"
)
data = rabpro.basin_stats.fetch_gee(urls, ["soc5"], ["system:index"])
res = gpd.GeoDataFrame(pd.concat([data, gdf], axis=1))
res.head()
[7]:
| soc5_mean | da_km2 | geometry | |
|---|---|---|---|
| 0 | 2407.976065 | 440.266532 | POLYGON ((-69.49583 44.92500, -69.50417 44.92500, -69.50502 44.92832, -69.50833 44.92917, -69.54... |