From 4571d2a589b93d51670ab52a21a4ce13794328e4 Mon Sep 17 00:00:00 2001
From: Axel Schnitger <axel.schnitger@aei.mpg.de>
Date: Mon, 20 Aug 2018 11:40:33 +0200
Subject: [PATCH] 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.
---
.gitignore | 1 +
README.md | 60 +++++++++++++++--------
function_gfr/readGFC.m | 105 +++++++++++++++++++++++++++++++++++++++++
gfr_parallel.m | 17 ++-----
4 files changed, 151 insertions(+), 32 deletions(-)
create mode 100644 function_gfr/readGFC.m
diff --git a/.gitignore b/.gitignore
index d607ce2..bcd060f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -49,3 +49,4 @@ octave-workspace
# End of https://www.gitignore.io/api/matlab,linux
results/
+input_data/
diff --git a/README.md b/README.md
index 34c17b1..42d6862 100644
--- a/README.md
+++ b/README.md
@@ -1,25 +1,51 @@
# GRACETOOLS
-
## Overview
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
-* [Usage](#usage)
-* [Features](#features)
-* [Content](#content)
-* [Acknowledgments](#acknowledgments)
-* [Contributors](#contributors)
+- [Usage](#usage)
+- [Features](#features)
+- [Files](#files)
+- [Acknowledgments](#acknowledgments)
+- [Contributors](#contributors)
## 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://wolke7.aei.mpg.de/index.php/s/DSoNjiclxl0S1UW)
+- [Download dataset](https://seafile.projekt.uni-hannover.de/f/3ee8fef6c6cb485489b4/?dl=1)
+
+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
@@ -43,16 +69,10 @@ input_data
├── EGM96.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.
-- 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
+## Files
| m-file | Description |
| :------------------------------ | :-------------------------------------------------------------------------------------------------------------------- |
@@ -94,7 +114,7 @@ The main mfiles to run are: `gfr_parallel.m` and `continue_gfr_par_from_iteratio
## Contributors
-- Neda Darbeheshti, neda.darbeheshti@aei.mpg.de
-- Florian Wöske, florian.woeske@zarm.uni-bremen.de
-- Axel Schnitger, axel.schnitger@aei.mpg.de
+- [Neda Darbeheshti](mailto:neda.darbeheshti@aei.mpg.de)
+- [Florian Wöske](mailto:florian.woeske@zarm.uni-bremen.de)
+- [Axel Schnitger](mailto:axel.schnitger@aei.mpg.de)
diff --git a/function_gfr/readGFC.m b/function_gfr/readGFC.m
new file mode 100644
index 0000000..5a9b90d
--- /dev/null
+++ b/function_gfr/readGFC.m
@@ -0,0 +1,105 @@
+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
+
diff --git a/gfr_parallel.m b/gfr_parallel.m
index 31749aa..fe055ce 100644
--- a/gfr_parallel.m
+++ b/gfr_parallel.m
@@ -1,7 +1,7 @@
%% GRACE gravity field recovery (GFR)
%
-% GFR_PARALLEL is the main m-file preforming the gravity field estimation. It is
-% optimised for parallel computing using local and global parameters for
+% GFR_PARALLEL is the main m-file preforming the gravity field estimation. It
+% is optimised for parallel computing using local and global parameters for
% 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
% the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an
@@ -13,9 +13,6 @@
% Neda Darbeheshti, AEI Hannover, 2018-07.
%%
-clear all
-% close all;
-% clc;
format longg;
%% variable inputs
@@ -58,14 +55,14 @@ GM = 0.3986004415E+15;
ae = 0.6378136300E+07;
% reference for comparison
-ggm05s = importGravityField('GGM05S.gfc');
+ggm05s = readGFC('GGM05S.gfc');
ggm05s = ggm05s(1:n_back,:); % cut to desired length
ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
% d_vec_i = ggm05s(:,3:4)-gravityf(:,3:4);
% d_vec_i = [ggm05s(:,1:2) d_vec_i];
%Load a-priori/ reference gravity field
-egm96 = importGravityField('EGM96.gfc');
+egm96 = readGFC('EGM96.gfc');
egm96 = egm96(1:n_back,:); % cut to desired length
egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
@@ -260,7 +257,7 @@ for ItrNo = 0:ItrNo_max
d_vec = [field_vec(:,1:2) d_vec];
dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);
- fprintf('done with iteration %d \n', ItrNo)
+ fprintf('done with iteration %d \n', ItrNo)
end
@@ -305,7 +302,3 @@ for i = 1:ItrNo_max+1
semilogy(2:lmaxcs,dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
end
-
-
-
-
--
GitLab