{
"cells": [
{
"cell_type": "markdown",
"id": "48ff48ae-a1bb-497a-b293-73de7b8e08fb",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# Compute Terrestrial Water Storage anomalies from GRACE-FO data\n",
"*Author: R. Rietbroek Feb 2024 (r.rietbroek@utwente.nl)*\n",
"\n",
"Since 2002, GRACE and its follow-on mission GRACE-FO have been providing monthly fields of so-called Stokes coefficients. These spherical harmonic coefficients represent gravity changes which are predominantly linked to the redistribution of water. The combined effect of those mass changes are captured as the [Essential Climate Variable Terrestrial Water Storage (TWS)](https://gcos.wmo.int/en/essential-climate-variables/tws).\n",
"\n",
"Stokes coefficients from the GRACE solutions can be converted to an equivalent water height, but several processing steps are often needed to obtain good results. This notebook provides an example for applying the following steps:\n",
"\n",
"* Obtain time-anomalies by subtracting a static gravity field from the monthly solutions\n",
"* *Todo:* Substitute the less accurate degree 2 coefficientsi with alternatives from a Satellite Laser Ranging solution (optional)\n",
"* *Todo:* Add degree 1 variations, which are needed to resolve for the Earth's wobbly around it's center of mass\n",
"* Filter the coefficients to remove high degreee noise\n",
"* Map (spherical harmonic synthesis) the results to a geographical area\n",
"\n",
"The data used in this example is a sample extracted from a [geoslurp enabled] PostGreSQl database, the file links and their metadata are organized in a sql database, for easier matching (e.g. in time), but the use of a database or sql file is not strictly necessary.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "c8dae3c8-e54b-4133-aca7-1261a471ca6e",
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"#Optionally enable autoreloading for development purposes. Note that this does not automagically reload the binary extensions\n",
"%load_ext autoreload\n",
"%autoreload 2\n",
"#also supress some warnings from pandas\n",
"import warnings\n",
"warnings.simplefilter(action='ignore', category=FutureWarning)"
]
},
{
"cell_type": "markdown",
"id": "ee109f2d-ea80-45d6-8052-7c8eacf3be63",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## 1. Load the necessary python modules"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "203afee4-b4e8-4a07-9b78-90a01160aef2",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"import xarray as xr\n",
"import shxarray\n",
"import numpy as np\n",
"import os\n",
"import gzip\n",
"import sqlite3\n",
"import tarfile\n",
"import matplotlib.pyplot as mpl"
]
},
{
"cell_type": "markdown",
"id": "c3921e86-e045-471d-b90c-a77cd3b401fd",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## 2. Load the sample data from the tarball (making use of the tables in the sqlite file)\n",
"shxarray can open files in various formats (e.g. [icgem](http://icgem.gfz-potsdam.de/home) or the common GSM-V6 format used in many GRACE processing centers). When installing shxarray, these methods become available to Xarray. The methods are exposed through choosing an appropriate engine in Xarray `open_dataset` function, e.g.:\n",
"\n",
"```\n",
"dsicgem=xr.open_dataset(icgemfilename,engine=\"icgem\")\n",
"dsgsmv6=xr.open_dataset(GSMV6filename,engine=\"gsmv6\")\n",
"\n",
"```\n",
"Gzipped compressed file are also allowed (ending in `.gz`).\n",
"\n",
"Alternatively, one can pass an open file descriptor to the open_dataset call. This latter way is used below."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "62b24000-8102-4c15-aadc-3b4d70772b35",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"# Helper class which aids in extracing the gzipped data files from the Tar file without the need to extract the data to the disk\n",
"class TarExtracter:\n",
" def __init__(self,tararchive):\n",
" \"\"\" Opens a tararchive and keep its open acces point available to functions within this class\"\"\"\n",
" self.tarar=tarfile.open(os.path.join(datadir,tararchive),'r:gz')\n",
" def getmember(self,member):\n",
" \"\"\"Extract a specific (gzipped) member from an open archive, and return an open file object\"\"\"\n",
" tarfileobj=self.tarar.extractfile(member)\n",
" return gzip.open(tarfileobj)\n",
" def close(self):\n",
" self.tarar.close()\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "2b97b1e3-286f-4dad-972c-744baa086ba6",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Tables present in the sqlite3 file ../../../sample-data/GRACEDataSample_2020.sql:\n",
"\tTable gracel2\n",
"\tTable static\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020001-2020031_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020001-2020031_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020032-2020060_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020032-2020060_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020061-2020091_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020061-2020091_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020092-2020121_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020092-2020121_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020122-2020152_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020122-2020152_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020153-2020182_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020153-2020182_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020183-2020213_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020183-2020213_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020214-2020244_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020214-2020244_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020245-2020274_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020245-2020274_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020275-2020305_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020275-2020305_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020306-2020335_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020306-2020335_GRFO_JPLEM_BC01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GSM-2_2020336-2020366_GRFO_JPLEM_BA01_0600.gz\n",
"Reading GRACEDataSample_2020_files.tgz:gracefol2_jpl_rl06/GAA-2_2020336-2020366_GRFO_JPLEM_BC01_0600.gz\n"
]
}
],
"source": [
"#The SQLITE database holds an inventory of the sample data\n",
"datadir='../../../sample-data'\n",
"sqlfile=os.path.join(datadir,'GRACEDataSample_2020.sql')\n",
"\n",
"conn=sqlite3.connect(sqlfile)\n",
"\n",
"print(f\"Tables present in the sqlite3 file {sqlfile}:\")\n",
"for res in conn.execute(\"SELECT name FROM sqlite_master WHERE type='table';\").fetchall():\n",
" print(f\"\\tTable {res[0]}\")\n",
"\n",
"tarname,staticgravmem=conn.execute(\"SELECT * FROM static LIMIT 1\").fetchone()[13].split(\":\")\n",
"tarar=TarExtracter(tarname)\n",
"\n",
"# read the time-invariable Stokes Coefficients\n",
"with tarar.getmember(staticgravmem) as fid:\n",
" dsstatic=xr.open_dataset(fid,engine=\"icgem\")\n",
"\n",
"# read the time variable Stokes coefficients\n",
"gsm=[]\n",
"gaa=[]\n",
"for res in conn.execute(f\"SELECT gsm,gaa FROM gracel2\").fetchall():\n",
" \n",
" print(f\"Reading {res[0]}\")\n",
" with tarar.getmember(res[0].split(\":\")[1]) as fid:\n",
" gsm.append(xr.open_dataset(fid,engine=\"gsmv6\"))\n",
" \n",
" print(f\"Reading {res[1]}\")\n",
" with tarar.getmember(res[1].split(\":\")[1]) as fid:\n",
" gaa.append(xr.open_dataset(fid,engine=\"gsmv6\"))\n",
" \n",
"dsgsm=xr.concat(gsm,dim=\"time\")\n",
"dsgaa=xr.concat(gaa,dim=\"time\")\n",
"\n",
"tarar.close()"
]
},
{
"cell_type": "markdown",
"id": "cb223110-a46d-4b08-bec3-b241bce4feae",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## 3. Create anomalies by subtracting the static gravity field from the monthly solutions"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "04a936c0-aa09-4d7f-9b0f-72d6f752bef3",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"dsgsm[\"dcnm\"]=dsgsm.cnm-dsstatic.cnm\n",
"\n",
"# Optional: restore (add) gaa atmospheric monthly averaged dealiasing product\n",
"# dsgsm[\"dcnm\"]=dsgsm.dcnm+dsgaa.cnm\n",
"nmax=dsgsm.sh.nmax"
]
},
{
"cell_type": "markdown",
"id": "f8bf6529-cea4-46b6-b4b7-d7ecc0aef70e",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## 4. Convert to Terrestrial Water Storage"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "ef61c7a9-5e70-4166-adc7-8978f1115712",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"shxarray-INFO: /home/roelof/.cache/shxarray_storage/Love/geoslurp_dump_llove.sql already exists, no need to download)\n"
]
},
{
"data": {
"text/html": [
"