GRACETOOLS

Overview

GRACETOOLS is a collection of MATLAB scripts that can be used for gravity field recovery using GRACE type satellite observations. There are two releases so far.

Features

The code is written using a batch least squares algorithm.
Three different fixed step numerical integration schemes are provided.
MATLAB parallel for-Loops (parfor) is used for calculation over arcs (days here).
This code can be easily modified to run on the user’s local parallel computations clusters.

GRACETOOLS is open source software and licensed under GPLv3. Please reference this article if you use GRACETOOLS:

Darbeheshti N., Wöske F., Weigelt M., McCullough C., Wu H. (2018) GRACETOOLS - GRACE gravity field recovery tools, Geosciences.

Requirements

MATLAB and the Parallel Computing Toolbox.

Download

Get the Scripts

Download GRACETOOLS as a zip file and extract it somewhere on your computer.

gracetools.zip

Or clone the repository with git clone:

$ git clone git@gitlab.aei.uni-hannover.de:geoq/gracetools.git

Get the Data

You need an example dataset to run these scripts. If you’re on Linux or MacOS use the get_input_data.sh script to get the dataset:

$ ./get_input_data.sh

If you prefer the manual way download the dataset:

Download dataset

Move the zip file into the input_data directory, extract it and delete the now unnecessary input_data.zip file. For completeness here are the terminal commands:

$ mv ~/Downloads/input_data.zip path/to/GRACETOOLS/input_data/
$ cd path/to/GRACETOOLS/input_data/
$ unzip input_data.zip
$ rm -f input_data.zip

The directory structure of input_data/ should look like this:

input_data
├── EGM96.gfc
├── GGM05S.gfc
├── MockData_do10_nod4
│   ├── 2005-05-01
│   │   ├── GNV1B_2005-05-01_A_02.asc
│   │   ├── GNV1B_2005-05-01_B_02.asc
│   │   └── KBR1B_2005-05-01_X_02.asc
│   ├── 2005-05-02
│   │   ├── GNV1B_2005-05-02_A_02.asc
│   │   ├── GNV1B_2005-05-02_B_02.asc
│   │   └── KBR1B_2005-05-02_X_02.asc
│   ├── 2005-05-03
│   │   ├── GNV1B_2005-05-03_A_02.asc
│   │   ├── GNV1B_2005-05-03_B_02.asc
│   │   └── KBR1B_2005-05-03_X_02.asc
│   └── 2005-05-04
│       ├── GNV1B_2005-05-04_A_02.asc
│       ├── GNV1B_2005-05-04_B_02.asc
│       └── KBR1B_2005-05-04_X_02.asc
└── MockData_do20_nod12
    ├── 2005-05-01
    │   ├── GNV1B_2005-05-01_A_02.asc
    │   ├── GNV1B_2005-05-01_B_02.asc
    │   └── KBR1B_2005-05-01_X_02.asc
    ├── 2005-05-02
    │   ├── GNV1B_2005-05-02_A_02.asc
    │   ├── GNV1B_2005-05-02_B_02.asc
    │   └── KBR1B_2005-05-02_X_02.asc
    ├── 2005-05-03
    │   ├── GNV1B_2005-05-03_A_02.asc
    │   ├── GNV1B_2005-05-03_B_02.asc
    │   └── KBR1B_2005-05-03_X_02.asc
    ├── 2005-05-04
    │   ├── GNV1B_2005-05-04_A_02.asc
    │   ├── GNV1B_2005-05-04_B_02.asc
    │   └── KBR1B_2005-05-04_X_02.asc
    ├── 2005-05-05
    │   ├── GNV1B_2005-05-05_A_02.asc
    │   ├── GNV1B_2005-05-05_B_02.asc
    │   └── KBR1B_2005-05-05_X_02.asc
    ├── 2005-05-06
    │   ├── GNV1B_2005-05-06_A_02.asc
    │   ├── GNV1B_2005-05-06_B_02.asc
    │   └── KBR1B_2005-05-06_X_02.asc
    ├── 2005-05-07
    │   ├── GNV1B_2005-05-07_A_02.asc
    │   ├── GNV1B_2005-05-07_B_02.asc
    │   └── KBR1B_2005-05-07_X_02.asc
    ├── 2005-05-08
    │   ├── GNV1B_2005-05-08_A_02.asc
    │   ├── GNV1B_2005-05-08_B_02.asc
    │   └── KBR1B_2005-05-08_X_02.asc
    ├── 2005-05-09
    │   ├── GNV1B_2005-05-09_A_02.asc
    │   ├── GNV1B_2005-05-09_B_02.asc
    │   └── KBR1B_2005-05-09_X_02.asc
    ├── 2005-05-10
    │   ├── GNV1B_2005-05-10_A_02.asc
    │   ├── GNV1B_2005-05-10_B_02.asc
    │   └── KBR1B_2005-05-10_X_02.asc
    ├── 2005-05-11
    │   ├── GNV1B_2005-05-11_A_02.asc
    │   ├── GNV1B_2005-05-11_B_02.asc
    │   └── KBR1B_2005-05-11_X_02.asc
    ├── 2005-05-12
    │   ├── GNV1B_2005-05-12_A_02.asc
    │   ├── GNV1B_2005-05-12_B_02.asc
    │   └── KBR1B_2005-05-12_X_02.asc
    └── states0_perfect.mat

Usage

Release 1

The main m-files in the first release (release1/) are:

gfr_parallel.m
save_vars2continue_itr.m
continue_gfr_par_from_iteration.m

First run gfr_parallel.m. Use save_vars2continue_itr.m to save all necessary variables in the results/ subdirectory. You can continue the estimation from the last iteration with continue_gfr_par_from_iteration.m.

Performance

With an older Intel i5 (4 cores) a test case with degree and order 10 with 4 days of observation, each iteration takes about 10 min.

Reference

Darbeheshti N., Wöske F., Weigelt M., McCullough C., Wu H. (2018) GRACETOOLS - GRACE gravity field recovery tools, Geosciences 2018, 8(9), 350; doi.org/10.3390/geosciences8090350.

Release 2

The main m-files in the second release (release2/) are:

gfr_eba.m - GFR with energy balance approach
gfr_pos.m - GFR with variational equations using positions
gfr_pos_rrho.m - GFR with variational equations using positions and range rates
gfr_rrho.m - GFR with variational equations using range rates
gfr_rrho_tickonov.m- GFR with variational equations using range rates and Tickonov regularization

Navigate to the release2/ subdirectory and run the scripts. Results will be stored in the results subdirectory.

Performance

With an older Intel i5 (4 cores) a test case with degree and order 20 with 12 days of observation, each iteration takes about 20 min.

Reference

Darbeheshti N., Wöske F., Weigelt M., Wu H., Mccullough C. (2020) Comparison of Spacewise and Timewise Methods for GRACE Gravity Field Recovery. In: Montillet JP., Bos M. (eds) Geodetic Time Series Analysis in Earth Sciences. Springer Geophysics. Springer, Cham, doi.org/10.1007/978-3-030-21718-1_10
IUGG Talk Montreal 2019

Files

From `function_gfr`

batch_processor_partitioned.m - batch processing algorithm for GRACE range-rate observations
cs2sc.m- converts the square containing spherical harmonics coefficients storage format into a rectangular format
cs2vec- rearranges a field of spherical harmonic coefficients in cs-or sc-format to a vector shape
deriv.m- calculates equation of motion and all partials
dv_geoidn.m- reads an Earth’s gravity field model and plots the degree variances
dv_geoidn_no_plot- reads an Earth’s gravity field model and gives the degree variances without plotting
grtenpshs.m- calculates the gradient and the tensor of the gravity field
hmat.m- makes the Hi_tilda matrix
importGravityField.m- import numeric data from a text file as a matrix
initplm.m- initializes a plm calculation
inter_sat_dist.m- calculates range and range rate from position and velocity
llpartialgradV.m- determines the partial derivative of the gradient of the gravity field w.r.t. the coefficients
multmatvek.m- mulitplicates a matrix and a vector
odeint.m- integrates one step with step size dt from t0 and initial values y0
odeint_abm.m- integrates one step with step size h from t0 and initial values y0
odeint_dp8.m- integrates one step with step size h from t0 and initial values y0 with Runge-Kutta (RK) integrator
odeint_dp8_abm_init.m- integrates one step with step size h from t0 and initial values y0 with Runge-Kutta (RK) integrator
odeint_rk4.m- integrates one step with step size h from t0 and initial values y0 with simple 4th order Runge-Kutta (RK) integrator
odeint_rk4_abm_init.m- integrates one step with step size h from t0 and initial values y0 with simple 4th order Runge-Kutta (RK) integrator
plm.m- fully normalized associated Legendre functions for a selected order M
plmsp.m- fully normalized associated Legendre functions for all degree and order but in a single points
readACC.m- Read GRACE ACC1B
readGFC.m- Read Gravity Field Coefficients (GFC file)
readGNV.m- Read GRACE GNV1B
readKBR.m- Read GRACE KBR1B
readSCA.m- Read GRACE SCA1B
Ri2e.m- returns the rotation matrix (only Z-rotation) from an inertial frame to an Earth-fixed coordinate system
Ri2e_dot.m- returns the derivative of the rotation matrix (only Z-rotation) from inertial frame to an Earth-fixed coordinate system
save_vars2continue_itr.m- save all important variables of a gravity estimation
vec2cs.m- rearranges a vector shaped set of spherical harmonic coefficients into cs-format

From `function_eba`

calculatePnm.m - Calculation of normalized Legendre Functions using stable recursion formulas
designmatrixn.m - Calculation of normalized Legendre Functions using stable recursion formulas
potential_EBAs.m - Calculates the potential along the GRACE simulated orbit
sortCoefficients.m - Sorts spherical harmonic coefficients

Contributors

Neda Darbeheshti, Florian Wöske, Axel Schnitger