Commit 2385c6bf authored by Fabian Wachsmann's avatar Fabian Wachsmann
Browse files

Delete use-case_count_summer_days_era5_and_cmip6.ipynb

parent 3da824de
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Count Annual Summer Days for a particular Location"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this Use Case you will learn the following:\n",
"- How to access a data file from DKRZ's CMIP6 archive\n",
"- How to access a remote file from the [Copernicus Climate Data Store](https://cds.climate.copernicus.eu/)\n",
"- How to count the annual number of summer days for a particular location in both data sets\n",
"- How to visualize the results\n",
"\n",
"\\\n",
"You will be using:\n",
"- [Intake](https://github.com/intake/intake) for finding the data the data\n",
"- [Xarray](http://xarray.pydata.org/en/stable/) for loading and processing the data\n",
"- [hvPlot](https://hvplot.holoviz.org/index.html) for visualizing the data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get formating done automatically according to style `black`\n",
"#%load_ext lab_black"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load Packages"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import intake # a general interface for loading data\n",
"\n",
"# import requests\n",
"import folium # visulization tool\n",
"import xarray as xr # handling lablled multi-dimensional arrays\n",
"\n",
"# import matplotlib\n",
"# from IPython.display import display\n",
"from ipywidgets import widgets # allows the use of widgets in Jupyer Notebook\n",
"from geopy.geocoders import Nominatim\n",
"import numpy as np # fundamenrtal package for scientific computing\n",
"import pandas as pd # data analysis and manipulation tool\n",
"import hvplot.pandas # visualization tool\n",
"import cdsapi # Climate Data Store API"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Choose Location and Year\n",
"If ambiguous, the more likely location will be chosen\n",
"\n",
"<a id='selection'></a>"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"place = widgets.Text(description=\"Enter place:\")\n",
"display(place)\n",
"\n",
"x = range(1950, 2014)\n",
"year = widgets.Dropdown(options=x, description=\"Select year: \", disabled=False,)\n",
"display(year)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Check if your Selection is correct"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"geolocator = Nominatim(user_agent=\"any_agent\")\n",
"location = geolocator.geocode(place.value)\n",
"\n",
"m = folium.Map(location=[location.latitude, location.longitude])\n",
"tooltip = location.latitude, location.longitude\n",
"folium.Marker([location.latitude, location.longitude], tooltip=tooltip).add_to(m)\n",
"display(m)\n",
"\n",
"print(location.address)\n",
"print((location.latitude, location.longitude))\n",
"print(\"Chosen year: \" + str(year.value))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Download ERA5 Data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Convert year and coordinates into strings for API request\n",
"coordinate_string = str(location.latitude) +'/' +str(location.longitude) +'/' +str(location.latitude) +'/' +str(location.longitude)\n",
"print(\"Coordinates: \" +coordinate_string)\n",
"year_string = str(year.value)\n",
"print(\"Year: \" +year_string)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# retrieve data via Climate Data Store Client\n",
"\n",
"c = cdsapi.Client()\n",
"\n",
"c.retrieve(\n",
" 'reanalysis-era5-single-levels',\n",
" {\n",
" 'variable': '2m_temperature',\n",
" 'product_type': 'reanalysis',\n",
" 'year': year_string, # our chosen year\n",
" 'month':[ # all months of this year\n",
" '01','02','03'\n",
" '04','05','06',\n",
" '07','08','09',\n",
" '10','11','12'\n",
" ],\n",
" 'day':[ # all days of the month\n",
" '01','02','03',\n",
" '04','05','06',\n",
" '07','08','09',\n",
" '10','11','12',\n",
" '13','14','15',\n",
" '16','17','18',\n",
" '19','20','21',\n",
" '22','23','24',\n",
" '25','26','27',\n",
" '28','29','30',\n",
" '31'\n",
" ],\n",
" 'time':[ # all hours of the day\n",
" '00:00','01:00','02:00',\n",
" '03:00','04:00','05:00',\n",
" '06:00','07:00','08:00',\n",
" '09:00','10:00','11:00',\n",
" '12:00','13:00','14:00',\n",
" '15:00','16:00','17:00',\n",
" '18:00','19:00','20:00',\n",
" '21:00','22:00','23:00'\n",
" ],\n",
" 'area' : coordinate_string, # North/ West/ South/ East. Default: global\n",
" \"format\": \"netcdf\", # save in NetCDF format\n",
" }, \"era5_2m_temperature.nc\" # save as ...\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Open File and Calculate Daily Maximum"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Open and access 2m temperature\n",
"t_xr = xr.open_dataset(\"era5_2m_temperature.nc\")['t2m']\n",
"\n",
"# have a look what's in there\n",
"t_xr"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# alternative approach\n",
"\n",
"tdaymax_xr = t_xr.groupby('time.dayofyear').max()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Convert to pandas series and calculate daily maximum by using resample function\n",
"t_s = pd.Series(t_xr[:,0,0].values, index=t_xr['time'].values)\n",
"tdaymax_s = t_s.resample('D').max()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare Reanalysis Grid Cell with chosen Location"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# draw map\n",
"m = folium.Map(location=[location.latitude, location.longitude], zoom_start=8)\n",
"\n",
"tooltip = location.latitude, location.longitude\n",
"folium.Marker(\n",
" [location.latitude, location.longitude],\n",
" tooltip=tooltip,\n",
" popup=\"Location selected by You\",\n",
").add_to(m)\n",
"\n",
"tooltip = float(t_xr['latitude'].values), float(t_xr['longitude'].values)\n",
"folium.Marker(\n",
" [t_xr['latitude'], t_xr['longitude']],\n",
" tooltip=tooltip,\n",
" popup=\"Location selected by ERA5\",\n",
").add_to(m)\n",
"\n",
"\n",
"# Define rectangle\n",
"rect_lat1_era5 = t_xr['latitude'][0] -0.125\n",
"rect_lon1_era5 = t_xr['longitude'][0] -0.125\n",
"rect_lat2_era5 = t_xr['latitude'][0] +0.125\n",
"rect_lon2_era5 = t_xr['longitude'][0] +0.125\n",
"#\n",
"\n",
"folium.Rectangle(\n",
" bounds=[[rect_lat1_era5, rect_lon1_era5], [rect_lat2_era5, rect_lon2_era5]],\n",
" color=\"#999900\",\n",
" fill=True,\n",
" fill_color=\"#ffff00\",\n",
" fill_opacity=0.2,\n",
").add_to(m)\n",
"m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load Intake Catalogue\n",
"We load the catalog descriptor with intake:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"col_url = \"/work/ik1017/Catalogs/mistral-cmip6.json\"\n",
"\n",
"col = intake.open_esm_datastore(col_url)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see what is in the intake catalog. The underlying data base is given as a panda dataframe which we can access with 'col.df'. col.df.head() shows us the first rows of the table:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"col.df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Browse Intake Catalogue and open Variable"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"query = dict(\n",
" source_id=\"MPI-ESM1-2-HR\", # here we choose Max-Plack Institute's Earth Sytem Model in high resolution\n",
" variable_id=\"tasmax\", # temperature at surface, maximum\n",
" table_id=\"day\", # daily maximum\n",
" experiment_id=\"historical\", # historical Simulation, 1850-2014\n",
" member_id=\"r10i1p1f1\",\n",
")\n",
"cat = col.search(**query)\n",
"\n",
"# Show query\n",
"cat.df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"open selected year"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"query_result_df = cat.df\n",
"query_result_df[\"start_year\"] = query_result_df[\"time_range\"].str[0:4].astype(int)\n",
"query_result_df[\"end_year\"] = query_result_df[\"time_range\"].str[9:13].astype(int)\n",
"query_result_df[\"selection\"] = (year.value - query_result_df[\"start_year\"] >= 0) & (\n",
" year.value - query_result_df[\"end_year\"] <= 0\n",
")\n",
"query_result_df[\"path\"].where(query_result_df[\"selection\"] == True)\n",
"selected_path_index = query_result_df.loc[query_result_df[\"selection\"] == True][\n",
" \"path\"\n",
"].index[0]\n",
"selected_path = query_result_df[\"path\"][selected_path_index]\n",
"\n",
"# Load Data\n",
"ds_tasmax = xr.open_dataset(selected_path)\n",
"tasmax_xr = ds_tasmax[\"tasmax\"]\n",
"\n",
"time_start = str(year.value) + \"-01-01T12:00:00.000000000\"\n",
"time_end = str(year.value) + \"-12-31T12:00:00.000000000\"\n",
"tasmax_year_xr = tasmax_xr.loc[time_start:time_end, :, :]\n",
"\n",
"# Show path for selected year\n",
"selected_path"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Let's have a look at the file\n",
"\n",
"tasmax_year_xr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare Model Grid Cell with chosen Location"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Find nearest model coordinate\n",
"\n",
"# First, find the index of the grid point nearest a specific lat/lon.\n",
"abslat = np.abs(tasmax_year_xr[\"lat\"] - location.latitude)\n",
"abslon = np.abs(tasmax_year_xr[\"lon\"] - location.longitude)\n",
"c = np.maximum(abslon, abslat)\n",
"\n",
"([xloc], [yloc]) = np.where(c == np.min(c))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"float(tasmax_year_xr[\"lat\"][yloc].values)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Draw map again\n",
"\n",
"m = folium.Map(location=[location.latitude, location.longitude], zoom_start=8)\n",
"\n",
"tooltip = float(tasmax_year_xr[\"lat\"][yloc].values), float(tasmax_year_xr[\"lon\"][xloc].values)\n",
"folium.Marker(\n",
" [tasmax_year_xr[\"lat\"][yloc], tasmax_year_xr[\"lon\"][xloc]],\n",
" tooltip=tooltip,\n",
" popup=\"Model Grid Cell Center\",\n",
").add_to(m)\n",
"\n",
"tooltip = location.latitude, location.longitude\n",
"folium.Marker(\n",
" [location.latitude, location.longitude],\n",
" tooltip=tooltip,\n",
" popup=\"Location selected by You\",\n",
").add_to(m)\n",
"\n",
"tooltip = float(t_xr['latitude'].values), float(t_xr['longitude'].values)\n",
"folium.Marker(\n",
" [t_xr['latitude'], t_xr['longitude']],\n",
" tooltip=tooltip,\n",
" popup=\"Location selected by ERA5\",\n",
").add_to(m)\n",
"\n",
"\n",
"# Define era5 rectangle\n",
"rect_lat1_era5 = t_xr['latitude'][0] -0.125\n",
"rect_lon1_era5 = t_xr['longitude'][0] -0.125\n",
"rect_lat2_era5 = t_xr['latitude'][0] +0.125\n",
"rect_lon2_era5 = t_xr['longitude'][0] +0.125\n",
"#\n",
"\n",
"folium.Rectangle(\n",
" bounds=[[rect_lat1_era5, rect_lon1_era5], [rect_lat2_era5, rect_lon2_era5]],\n",
" color=\"#999900\",\n",
" fill=True,\n",
" fill_color=\"#ffff00\",\n",
" fill_opacity=0.2,\n",
").add_to(m)\n",
"\n",
"# Define model rectangle\n",
"rect_lat1_model = (tasmax_year_xr[\"lat\"][yloc - 1] + tasmax_year_xr[\"lat\"][yloc]) / 2\n",
"rect_lon1_model = (tasmax_year_xr[\"lon\"][xloc - 1] + tasmax_year_xr[\"lon\"][xloc]) / 2\n",
"rect_lat2_model = (tasmax_year_xr[\"lat\"][yloc + 1] + tasmax_year_xr[\"lat\"][yloc]) / 2\n",
"rect_lon2_model = (tasmax_year_xr[\"lon\"][xloc + 1] + tasmax_year_xr[\"lon\"][xloc]) / 2\n",
"#\n",
"\n",
"folium.Rectangle(\n",
" bounds=[[rect_lat1_model, rect_lon1_model], [rect_lat2_model, rect_lon2_model]],\n",
" color=\"#009900\",\n",
" fill=True,\n",
" fill_color=\"#00ff00\",\n",
" fill_opacity=0.2,\n",
").add_to(m)\n",
"\n",
"\n",
"m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Draw Temperature Time Series and Count Summer days"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tdaymax_s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tasmax_year_place_xr = tasmax_year_xr[:, yloc, xloc] - 273.15\n",
"tasmax_year_place_df = pd.DataFrame(index = tasmax_year_place_xr['time'].values, columns = ['Temperature', 'Summer Day Threshold'])\n",
"#\n",
"tasmax_year_place_df.loc[:, 'Model Temperature'] = tasmax_year_place_xr.values\n",
"tasmax_year_place_df.loc[:, 'Reanalysis Temperature'] = tdaymax_s.values -273.15\n",
"tasmax_year_place_df.loc[:, 'Summer Day Threshold'] = 25\n",
"\n",
"tasmax_year_place_df.hvplot.line(y=['Model Temperature', 'Reanalysis Temperature', 'Summer Day Threshold'], \n",
" value_label='Temperature in °C', legend='bottom', title='Daily maximum Temperature near Surface for ' +place.value, height=500, width=620)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"no_summer_days_model = tasmax_year_place_xr[tasmax_year_place_xr > 25].size\n",
"no_summer_days_reanalysis = tdaymax_s[tdaymax_s -273.15 > 25].size\n",
"print(\"Model: In \" +str(year.value) +\" \" +place.value +\" had \" + str(no_summer_days_model)+\" summer days.\")\n",
"print(\"Reanalysis: In \" +str(year.value) +\" \" +place.value +\" had \" + str(no_summer_days_reanalysis)+\" summer days.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Try another location and year](#selection)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "test_env",
"language": "python",
"name": "test_env"
},
"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.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment