Skip to content
Snippets Groups Projects
Commit 4571d2a5 authored by Axel Schnitger's avatar Axel Schnitger
Browse files

Several modifications and fixes.

- Add readGFC function to be able to read gfc files with headers.
- Change link to the example dataset and add a description.
- Brush up README.md
- Update ignore list.
parent 677a5ab1
No related branches found
No related tags found
No related merge requests found
...@@ -49,3 +49,4 @@ octave-workspace ...@@ -49,3 +49,4 @@ octave-workspace
# End of https://www.gitignore.io/api/matlab,linux # End of https://www.gitignore.io/api/matlab,linux
results/ results/
input_data/
# GRACETOOLS # GRACETOOLS
## Overview ## Overview
GRACETOOLS is an archive of MATLAB software that can be used for gravity field recovery using GRACE type satellite observations. GRACETOOLS is an archive of MATLAB software that can be used for gravity field recovery using GRACE type satellite observations.
**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.
## Content ## Content
* [Usage](#usage) - [Usage](#usage)
* [Features](#features) - [Features](#features)
* [Content](#content) - [Files](#files)
* [Acknowledgments](#acknowledgments) - [Acknowledgments](#acknowledgments)
* [Contributors](#contributors) - [Contributors](#contributors)
## Usage ## Usage
To be able to run the program download the following dataset:
To be able to run the program download the following dataset and extract it into the `input_data` directory: - [Download dataset](https://seafile.projekt.uni-hannover.de/f/3ee8fef6c6cb485489b4/?dl=1)
- [Download dataset](https://wolke7.aei.mpg.de/index.php/s/DSoNjiclxl0S1UW)
Move the dataset into the `input_data` directory and extract it. Then delete
the `input_data.zip` file. In a terminal the commands look like this from the
GRACETOOLS main directory:
```
$ mv ~/Downloads/input_data.zip input_data/
$ cd input_data/
$ unzip input_data.zip
$ rm -f input_data.zip
```
The directory structure of `input_data` should look lie this: Again from the GRACETOOLS main directory you can check the directory structure
of the input_data/ directory with the following command:
```
$ tree input_data/
```
The output should look like this:
``` ```
input_data input_data
...@@ -43,16 +69,10 @@ input_data ...@@ -43,16 +69,10 @@ input_data
├── EGM96.gfc ├── EGM96.gfc
└── GGM05S.gfc └── GGM05S.gfc
``` ```
The main mfiles to run are: `gfr_parallel.m` and `continue_gfr_par_from_iteration.m`.
## Features The main m-files to run are: `gfr_parallel.m` and `continue_gfr_par_from_iteration.m`.
- The code is written using a batch least squares algorithm. ## Files
- 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.
## Content
| m-file | Description | | m-file | Description |
| :------------------------------ | :-------------------------------------------------------------------------------------------------------------------- | | :------------------------------ | :-------------------------------------------------------------------------------------------------------------------- |
...@@ -94,7 +114,7 @@ The main mfiles to run are: `gfr_parallel.m` and `continue_gfr_par_from_iteratio ...@@ -94,7 +114,7 @@ The main mfiles to run are: `gfr_parallel.m` and `continue_gfr_par_from_iteratio
## Contributors ## Contributors
- Neda Darbeheshti, neda.darbeheshti@aei.mpg.de - [Neda Darbeheshti](mailto:neda.darbeheshti@aei.mpg.de)
- Florian Wöske, florian.woeske@zarm.uni-bremen.de - [Florian Wöske](mailto:florian.woeske@zarm.uni-bremen.de)
- Axel Schnitger, axel.schnitger@aei.mpg.de - [Axel Schnitger](mailto:axel.schnitger@aei.mpg.de)
function gfc = readGFC(filename, verbosity)
%------------------------------------------------------------------------------
%
% READGFC(FILENAME, VERBOSITY) loads a global gravity field model in form of a
% spherical harmonic expansion (gcf-file).
%
% Input:
%
% filename: path to a .gfc file
% verbosity: if set to true additional output to command window. Default is
% false.
%
% Output:
%
% gfc: matrix containing l,m,C,S
%
% Usage:
%
% readGFC('filename')
% readGFC('filename', true)
%
%------------------------------------------------------------------------------
% Authors:
% Axel Schnitger, AEI, Leibniz University Hannover
%------------------------------------------------------------------------------
% Revision History
% 2018-08-20: AS, intitial version
%------------------------------------------------------------------------------
% Input arguments
if nargin < 2 || verbosity == false
verbosity = false;
elseif verbosity
tic
fprintf('\nStart reading:\t%s\n', filename)
end
fileID = fopen(filename);
line = fgetl(fileID);
if contains(line, 'gfc')
hasHeader = false;
else
hasHeader = true;
end
formatSpec = '%s %f %f %f %f %f %f';
if hasHeader
headerLines = 1;
while ~contains(line, 'end_of_head')
line = fgetl(fileID);
if isempty(line)
headerLines = headerLines + 1;
continue
end
keyword = textscan(line,'%s',1);
if(strcmp(keyword{1}, 'max_degree'))
tmp = textscan(line,'%s%d',1);
maxDegree = tmp{2};
end
if(strcmp(keyword{1}, 'radius'))
tmp = textscan(line,'%s%f',1);
radius = tmp{2};
end
if(strcmp(keyword{1}, 'earth_gravity_constant'))
tmp = textscan(line,'%s%f',1);
GM = tmp{2};
end
headerLines = headerLines + 1;
end
fclose(fileID);
fileID = fopen(filename);
gfcData = textscan(fileID, formatSpec, 'HeaderLines', headerLines);
fclose(fileID);
elseif ~hasHeader
headerLines = 0;
fileID = fopen(filename);
gfcData = textscan(fileID, formatSpec);
fclose(fileID);
maxDegree=max(gfcData{2});
end
if verbosity
if hasHeader
fprintf('Headerlines:\t%u\n', headerLines)
fprintf('Size data:\t%u\n', length(gfcData{3}))
fprintf('max. Degree:\t%u\n', maxDegree)
fprintf('Radius:\t\t%f\n', radius)
fprintf('GM:\t\t%f\n', GM)
elseif ~hasHeader
fprintf('No header found.\n');
fprintf('Size data:\t%u\n', length(gfcData{3}))
fprintf('max. Degree:\t%u\n', maxDegree)
end
toc
end
gfc = [gfcData{2}, gfcData{3}, gfcData{4}, gfcData{5}];
end
%% GRACE gravity field recovery (GFR) %% GRACE gravity field recovery (GFR)
% %
% GFR_PARALLEL is the main m-file preforming the gravity field estimation. It is % GFR_PARALLEL is the main m-file preforming the gravity field estimation. It
% optimised for parallel computing using local and global parameters for % is optimised for parallel computing using local and global parameters for
% different arcs. All necessary parameters and data are specified in the % different arcs. All necessary parameters and data are specified in the
% following. You have to set a number of iterations you want to perform. With % following. You have to set a number of iterations you want to perform. With
% the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an % the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an
...@@ -13,9 +13,6 @@ ...@@ -13,9 +13,6 @@
% Neda Darbeheshti, AEI Hannover, 2018-07. % Neda Darbeheshti, AEI Hannover, 2018-07.
%% %%
clear all
% close all;
% clc;
format longg; format longg;
%% variable inputs %% variable inputs
...@@ -58,14 +55,14 @@ GM = 0.3986004415E+15; ...@@ -58,14 +55,14 @@ GM = 0.3986004415E+15;
ae = 0.6378136300E+07; ae = 0.6378136300E+07;
% reference for comparison % reference for comparison
ggm05s = importGravityField('GGM05S.gfc'); ggm05s = readGFC('GGM05S.gfc');
ggm05s = ggm05s(1:n_back,:); % cut to desired length ggm05s = ggm05s(1:n_back,:); % cut to desired length
ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4)); ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
% d_vec_i = ggm05s(:,3:4)-gravityf(:,3:4); % d_vec_i = ggm05s(:,3:4)-gravityf(:,3:4);
% d_vec_i = [ggm05s(:,1:2) d_vec_i]; % d_vec_i = [ggm05s(:,1:2) d_vec_i];
%Load a-priori/ reference gravity field %Load a-priori/ reference gravity field
egm96 = importGravityField('EGM96.gfc'); egm96 = readGFC('EGM96.gfc');
egm96 = egm96(1:n_back,:); % cut to desired length egm96 = egm96(1:n_back,:); % cut to desired length
egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4)); egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
...@@ -305,7 +302,3 @@ for i = 1:ItrNo_max+1 ...@@ -305,7 +302,3 @@ for i = 1:ItrNo_max+1
semilogy(2:lmaxcs,dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str); semilogy(2:lmaxcs,dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
end end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment