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.

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:

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