-
Notifications
You must be signed in to change notification settings - Fork 1
/
RunDemo.m
159 lines (147 loc) · 9.59 KB
/
RunDemo.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
% Gravimetric quasigeoid estimation
%
% RunDemo is the main m-file preforming the gravimetric quasiGeoid estimation.
%
% Main parameters to be defined:
% - Grid/Tiling Parameters
% - DEM (Digital Elevation Model) data
% - Gravity data
% - Gravity Gradiometry data
% - Covariance function parameters
% - Topo condensation parameters
% - GGM (Global Gravity Model) reference signal
% - Coastline data
% - Levelling data comparisons
% - Output directories
%
% Main functions:
% - ImportAndFormatData: import all different data sets needed for LSC.
% - computeTerrainEffect: calculates the Residual Terrain Model.
% - computeGravimetryLSC: runs the LSC (Least Squares Collocation) in blocks.
% - MosaicTiles: mosaic the tiles with weights and save final result in tif
%
% Written by Jack McCubbine
% Last updated by Neda Darbeheshti
% Geoscience Australia, 2024-03.
%%
close all
clear
% Turn off the irritating warnings.
warning off
%% Add the path to the function files.
addpath('functions');
%% Grid/Tiling Parameters
% Tiling Parameters - fixed for each computation
GRID_PARA.buffer=1;% degs. The x/y extent to extract data around the tile centre. .75
GRID_PARA.buffer2=0.5;% degs. The x/y tile extents that are kept - where the good data are.
GRID_PARA.STEP=0.5;% The step size. This must be less than buffer2 to avoid gaps in the final grid.
GRID_PARA.filterSize=15;% filter size for spatial grid weight, this value is from experiment for tiles of one degree
GRID_PARA.filterRadius=10; % filter radius for spatial grid weight, this value is from experiment for tiles of one degree
% Grid extents - ensure these values are in GRID_PARA.STEP degree value increments.
% Boundary for computation
GRID_PARA.MINLONG=141.5;
GRID_PARA.MAXLONG=142;
GRID_PARA.MINLAT=-32.5;
GRID_PARA.MAXLAT=-32;
%% DEM data - N.B. the dem is used to specify the grid nodes.
DEM_PARA.filename='Data\DEM\AUSDEM1min.xyz';
DEM_PARA.num_cols=4861; % original size, to be changed, see below
DEM_PARA.num_rows=3181; % original size, to be changed, see below
%% Gravity data
% this collates all of the gravity and position data into one matlab array.
GRAV_PARA.filename='Data\processedData\Terrestrial_Gravity.mat';
GRAV_PARA.filename1=[];%'Data\GRAVITY\Xcalibur_Gravity.mat';
GRAV_PARA.TypeB=1;% This is a Type B uncertainty value (in mGal) which is added to the uncertainty values.
GRAV_PARA.Grav_Faye_TypeB=3;
%% Gravity Gradiometry data
GRAV_GRAD_PARA.filename='Data\GRAVITY_GRAD\Xcalibur_FVD_GDD.mat';
GRAV_GRAD_PARA.TypeB=10^(-5);% This is a Type B uncertainty value (in mGal/m) which is added to the uncertainty values.
GRAV_GRAD_PARA.avail=false;
%% Covariance function parameters
COV_PARA.Compute_Empircal_COV_Dec=3; % Decimation factor for empirical covariance estimates. e.g. 1 is no decimation, 2 drops 50% of the data etc. see sph_empcov for logic.
COV_PARA.Fit_Empircal_COV='auto';%'auto';% process to fit covariance N & M function values 'man' for manual to fit them on the cmd line or 'auto' , '' to just use what you supply here.
COV_PARA.FitEmpiricalCOVNSearch=[21600,1,21600];% Start, step, stop parameter sweep values for N parameter - if auto
COV_PARA.FitEmpiricalCOVMSearch=[200,20,300];% Start, step, stop parameter sweep values for M parameter - if auto
COV_PARA.N=10800;% max Legender polynonial of cov func.
COV_PARA.M=200;% min Legender polynonial of cov func.
COV_PARA.width=3;% Size of precomputed cov function in degrees - must be larger the the distance between any two points on a tile.
COV_PARA.res=30/3600; % Resolution of the covariance function
COV_PARA.COV_COMPUTED_Tilewise=true;% This recomputes the covariance function for each tile.
COV_PARA.Airbornedataonly=false;% Only use airborne data in establishing Covariance parameters - good to use if we are using EGM2008 as the references as terrestrial data are not independent.
COV_PARA.COVPlot=true;% true plots progress, false turns this off.
%% Topo condensation parameters
Topo_PARA.Corr=true;% MAKE SURE YOU TURN THIS ON!!!
Topo_PARA.TopoPlot=true;% true plots progress, false turns this off.
Topo_PARA.Density=2.67;% Assumed density in g/cm^3.
Topo_PARA.Depth=0;% Condensation layer depth. 0 is on the geoid
Topo_PARA.Rad=1;% Radius out to which to compute the effects in degress.
Topo_PARA.RTM=[50,10,300];%[1000,10,2160]for egm%[0,10,300]% Range of SHM degree filter parameters (min, step, max) explored when running RTM calculations.
%about 1080 for EGM,we need to fix the filter by factor of 2
%% GGM reference signal
% First run e.g. ./Data/GGM/RunIsGrafLab_Topo_Surf_EGM2008.m
GGM_PARA.filename='Data\processedData\GOCE_N300_For_Gridded_Int24_2height.mat';
%% Coastline data
COAST_PARA.filename='Data\COASTLINE\CoastAus.mat';
%% Levelling data comparisons
LEVELLING_PARA.Lev_eval=true;% If true, the levelling data are compared to the geoid as its computed.
LEVELLING_PARA.filename='Data\GPS_LEVELLING\Lev_NSW_NG.mat';%'Data/GPS_LEVELLING/Lev_CARS.mat';% The format of these data needs to be an array with rows [Long,Lat,h-H].
LEVELLING_PARA.Plot_Stats=true;% If true, the levelling data are compared to the geoid as its computed.
LEVELLING_PARA.Compare_To_Existing_Model=true;% If true, the levelling data are also compared to another existing geoid as its computed.
LEVELLING_PARA.Existing_Model='Data\EXISTING_GEOID_MODELS\AGQG20221120.mat';% File location of the existing model.
LEVELLING_PARA.max_diff=0.15;% Threshold for an outlier with the GNSS-levelling
%% Output
OUTPUT_PARA.Grids_name='outputs/Grids/';
OUTPUT_PARA.Tiles_dir_name='outputs/ResidualTiles/';
OUTPUT_PARA.PLOT_GRIDS=true;% A gridded solution is plotted and output as well as the tiles.
OUTPUT_PARA.plotsFolder=['outputs/plots/',date,'demo'];
% Keep the computer awake
keepawake=true;% Setting this to true wiggles the mouse every so often so the compute doesnt go to sleep.
%% Run the LSC code
% this is how demo data is created
% disp('1/4 ..........................importAndFormatData is running ')
% [Gravity6D,gravGradFiltered,DEM3D,ZDEM_griddedInterpolant,LongDEMmatrix,LatDEMmatrix,...
% GravityGGM_griddedInterpolant,ZetaGGM_griddedInterpolant,LevellingData3D,...
% ZetaRef_griddedInterpolant,GridRef3D,Coastline,DEM_PARA]=importAndFormatData...
% (GRID_PARA,DEM_PARA,GRAV_PARA,Topo_PARA,COAST_PARA,LEVELLING_PARA,GGM_PARA,GRAV_GRAD_PARA);
%
% save([OUTPUT_PARA.Grids_name,'oneTileData141-32','.mat'],...
% 'Gravity6D','DEM3D','LongDEMmatrix','LatDEMmatrix','LevellingData3D','GridRef3D','Coastline')
%
% save([OUTPUT_PARA.Grids_name,'oneTileGriddedInterpolant141-32','.mat'],...
% 'ZDEM_griddedInterpolant','GravityGGM_griddedInterpolant','ZetaGGM_griddedInterpolant','ZetaRef_griddedInterpolant')
disp('1/4 ..........................import demo data')
dataDemo=importdata('Data\oneTileData141-32.mat');
interpolantDemo=importdata('Data\oneTileGriddedInterpolant141-32.mat');
% Plot input data:
plotCustomScatter(dataDemo.DEM3D(:,1),dataDemo.DEM3D(:,2),dataDemo.DEM3D(:,3),GRID_PARA,'DEM','m',dataDemo.Coastline,OUTPUT_PARA.plotsFolder)
%caxis([-1 584])% always check plotCustomScatter for %caxis
plotCustomScatter(dataDemo.Gravity6D(:,1),dataDemo.Gravity6D(:,2),dataDemo.Gravity6D(:,3),GRID_PARA,'GravityHeight','m',dataDemo.Coastline,OUTPUT_PARA.plotsFolder)
%caxis([-1 584])
plotCustomScatter(dataDemo.Gravity6D(:,1),dataDemo.Gravity6D(:,2),dataDemo.Gravity6D(:,4),GRID_PARA,'Gravity','mGal',dataDemo.Coastline,OUTPUT_PARA.plotsFolder)
%caxis([-44 60])
disp('2/4 ..........................computeTerrainEffect is running')
[fullTopoCorrectedGravityPoint,longwaveTopo_griddedInterpolant,fullTopo_griddedInterpolant]=computeTerrainEffect(GRID_PARA, ...
Topo_PARA,dataDemo.Gravity6D,interpolantDemo.GravityGGM_griddedInterpolant,dataDemo.DEM3D,interpolantDemo.ZDEM_griddedInterpolant, ...
dataDemo.LongDEMmatrix,dataDemo.LatDEMmatrix,dataDemo.Coastline,OUTPUT_PARA.plotsFolder);
save([OUTPUT_PARA.Grids_name,'TerrainEffects','.mat'],'fullTopoCorrectedGravityPoint','longwaveTopo_griddedInterpolant','fullTopo_griddedInterpolant')
% if TerrainEffects have been saved, uncomment
% TE=importdata([OUTPUT_PARA.Grids_name,'TerrainEffects.mat']);
DEM_PARA.num_rows=size(dataDemo.LatDEMmatrix,1);
DEM_PARA.num_cols=size(dataDemo.LatDEMmatrix,2);
tic
disp('3/4 ..........................computeGravimetryLSC is running')
computeGravimetryLSC(GRID_PARA,COV_PARA,DEM_PARA,GRAV_PARA,OUTPUT_PARA,dataDemo.GridRef3D,fullTopoCorrectedGravityPoint, ...
interpolantDemo.GravityGGM_griddedInterpolant,interpolantDemo.ZDEM_griddedInterpolant,fullTopo_griddedInterpolant, ...
longwaveTopo_griddedInterpolant,Topo_PARA.Density)
toc
disp('4/4 ..........................mosaicTiles is running')
geomGravGeoidDiff = mosaicTiles(GRID_PARA,DEM_PARA,OUTPUT_PARA,dataDemo.LevellingData3D,dataDemo.LongDEMmatrix,dataDemo.LatDEMmatrix, ...
interpolantDemo.ZetaRef_griddedInterpolant,interpolantDemo.GravityGGM_griddedInterpolant,interpolantDemo.ZetaGGM_griddedInterpolant,dataDemo.Coastline);
% % test Parallel pull for tilewise LSC
% tic
%
% disp('3/4 ..........................computeParallelGravimetryLSC')
% computeParallelGravimetryLSC(GRID_PARA,COV_PARA,DEM_PARA,GRAV_PARA,OUTPUT_PARA,dataDemo.GridRef3D,fullTopoCorrectedGravityPoint, ...
% interpolantDemo.GravityGGM_griddedInterpolant,interpolantDemo.ZDEM_griddedInterpolant,fullTopo_griddedInterpolant, ...
% longwaveTopo_griddedInterpolant,Topo_PARA.Density)
% toc