Skip to content

Commit

Permalink
add fapar
Browse files Browse the repository at this point in the history
  • Loading branch information
@malmans2
malmans2 committed Jul 19, 2024
1 parent 2596914 commit dc59cff
Showing 1 changed file with 253 additions and 0 deletions.
253 changes: 253 additions & 0 deletions notebooks/wp5/satellite_fapar_completeness.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,253 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"# FAPAR satellite data completeness"
]
},
{
"cell_type": "markdown",
"id": "1",
"metadata": {},
"source": [
"## Import packages"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2",
"metadata": {},
"outputs": [],
"source": [
"import fsspec\n",
"import geopandas\n",
"import matplotlib.pyplot as plt\n",
"import shapely.geometry\n",
"from c3s_eqc_automatic_quality_control import download, utils\n",
"\n",
"plt.style.use(\"seaborn-v0_8-notebook\")"
]
},
{
"cell_type": "markdown",
"id": "3",
"metadata": {},
"source": [
"## Define variables"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4",
"metadata": {},
"outputs": [],
"source": [
"# Request parameters\n",
"year_start = 2014\n",
"year_stop = 2014\n",
"nominal_days = [3, 13, 21, 23, 24]\n",
"variables = [\"fapar\"]\n",
"\n",
"# Region\n",
"lon_slice = slice(-13, 35)\n",
"lat_slice = slice(72, 30)\n",
"\n",
"shapefile_url = \"https://figshare.com/ndownloader/files/23392280\""
]
},
{
"cell_type": "markdown",
"id": "5",
"metadata": {},
"source": [
"## Set the data request"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6",
"metadata": {},
"outputs": [],
"source": [
"# Define request\n",
"collection_id = \"satellite-lai-fapar\"\n",
"request = {\n",
" \"satellite\": \"proba\",\n",
" \"sensor\": \"vgt\",\n",
" \"horizontal_resolution\": \"1km\",\n",
" \"product_version\": \"V2\",\n",
" \"year\": [str(year) for year in range(year_start, year_stop + 1)],\n",
" \"month\": [f\"{month:02d}\" for month in range(1, 12 + 1)],\n",
" \"nominal_day\": [f\"{day:02d}\" for day in nominal_days],\n",
" \"format\": \"zip\",\n",
" \"variable\": variables,\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "7",
"metadata": {},
"source": [
"## Define function to compute missing values count"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8",
"metadata": {},
"outputs": [],
"source": [
"def compute_missing_values_count(ds):\n",
" ds.rio.set_spatial_dims(x_dim=\"longitude\", y_dim=\"latitude\", inplace=True)\n",
" ds.rio.write_crs(\"epsg:4326\", inplace=True)\n",
"\n",
" da_mvc = ds[\"fAPAR\"].isnull().sum(\"time\") / ds.sizes[\"time\"] * 100\n",
" da_mvc.attrs[\"long_name\"] = \"Missing values\"\n",
" da_mvc.attrs[\"units\"] = \"%\"\n",
"\n",
" return da_mvc.to_dataset(name=\"fAPAR\")"
]
},
{
"cell_type": "markdown",
"id": "9",
"metadata": {},
"source": [
"## Download and preprocess data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10",
"metadata": {},
"outputs": [],
"source": [
"# Download and cache\n",
"ds = download.download_and_transform(\n",
" collection_id,\n",
" request,\n",
" transform_func=utils.regionalise,\n",
" transform_func_kwargs={\n",
" \"lon_slice\": lon_slice,\n",
" \"lat_slice\": lat_slice,\n",
" },\n",
" chunks={\"year\": 1, \"nominal_day\": 1, \"variable\": 1},\n",
" cached_open_mfdataset_kwargs={\"combine\": \"nested\", \"concat_dim\": \"time\"},\n",
")\n",
"\n",
"# Shapefile\n",
"with fsspec.open(f\"simplecache::{shapefile_url}\") as file:\n",
" world_shape = geopandas.read_file(file, layer=\"Continents\")"
]
},
{
"cell_type": "markdown",
"id": "11",
"metadata": {},
"source": [
"## Define plotting function"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12",
"metadata": {},
"outputs": [],
"source": [
"def imshow_and_hist(da, shape):\n",
" \"\"\"Plot map and histogram side-by-side.\n",
"\n",
" Parameters\n",
" ----------\n",
" da: DataArray\n",
" DataArray to plot\n",
" shape: GeoDataFrame\n",
" Geopandas object with polygons\n",
"\n",
" Returns\n",
" -------\n",
" figure, axes\n",
" \"\"\"\n",
" fig, (ax_imshow, ax_hist) = plt.subplots(\n",
" 1, 2, figsize=[10, 5], gridspec_kw={\"width_ratios\": [3, 2]}\n",
" )\n",
"\n",
" da = da.rio.clip(\n",
" shape.geometry.apply(shapely.geometry.mapping),\n",
" shape.crs,\n",
" drop=True,\n",
" )\n",
" da.plot.imshow(ax=ax_imshow)\n",
" ax_imshow.set_title(\"Map\")\n",
"\n",
" da.plot.hist(bins=50, ax=ax_hist)\n",
" ax_hist.set_ylabel(\"Frequency\")\n",
" ax_hist.yaxis.set_label_position(\"right\")\n",
" ax_hist.yaxis.tick_right()\n",
"\n",
" # Compute and show no data percentage\n",
" missing_data_perc = (da == 100).sum() / da.notnull().sum() * 100\n",
" ax_hist.set_title(\n",
" f\"Percentage of pixels with no data respect the total: {float(missing_data_perc):f} %\"\n",
" )\n",
"\n",
" fig.suptitle(\", \".join(list(shape.CONTINENT)))\n",
" return fig, (ax_imshow, ax_hist)"
]
},
{
"cell_type": "markdown",
"id": "13",
"metadata": {},
"source": [
"## Plot Europe map"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "14",
"metadata": {},
"outputs": [],
"source": [
"da_mvc = ds[\"fAPAR\"].isnull().sum(\"time\") / ds.sizes[\"time\"] * 100\n",
"da_mvc.attrs[\"long_name\"] = \"Missing values\"\n",
"da_mvc.attrs[\"units\"] = \"%\"\n",
"\n",
"da_mvc.rio.write_crs(\"epsg:4326\", inplace=True)\n",
"imshow_and_hist(da_mvc, world_shape[world_shape.CONTINENT == \"Europe\"])"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}

0 comments on commit dc59cff

Please sign in to comment.