Geoid modelling ce678 project
- Jainam
- Vyush
- Sachin
This is a guide on how to use this software We use classical remove-restore method to find geoid height in this code.
If stuck anywhere use this guide or type help corresponding to the function in Matlab.
There is a also a file project.ipynb
that we have created to demonstrate our code on Nebraska,USA.
Also project.pdf
is the project report file.
-
Download the airborne gravity data from here by choosing any state tile.
-
Use the
import_airborne_gravity_data.m
function to import the data from the file of gravity data. -
Save the matrix into vectors . Ex -
T = import_airborne_gravity_data('file_path');`
latitude = T(:,1);
longitude = T(:,2);
elipsoidal_height = T(:,4);
gravity_filtered = T(:,3);
-
Download the geoid heights of GGM model from here
-
You can import them with import_geoid_heights
-
Code can be written as -
nggm = import_geoid_heights('filepath',latitude,longitude)
-
Calculate the Orthometric height by adding
nggm
toellisoidal_height
using this function-orthometric_height = calc_orthometric_height(elipsoidal_height,nggm)
Calculate the free air anomaly using the formula
Code -free_air_anomily = calc_free_air_anomaly(gravity_filtered,latitude,Ellipsoid,orthometric_height)
Here Ellipsoid
is a struct
used for storing important constants of Earth.For example ,You can use it by initiating as -
Ellipsoid = make_Ellipsoid('WGS84')
From the free-air gravity anomaly, remove the effect of a global gravity field model, a long- wavelength gravity anomaly
The resultant reduced gravity anomaly contains only the medium and the short wave- lengths.Download gravity anomaly of any GGM model from here.
TIPS on choosing-
- Choose latitude and longitude range within the tile of our gravity airborne data we got in 1st step.
- The website doesent allow resolution less than about 0.01. So make sure choose appropriate resolution so that code does not take long time to execute.
Code-
- To import the data
del_g_ggm = import_grav_anomaly('filepath',latitude,longitude)
del_g_s_and_mw = remove_g_ggm(free_air_anomily,del_g_ggm);
A correction to account for the gravitational attraction of the attraction must be applied
Code-
g_atm_and_s_mw = remove_g_atm(del_g_s_and_mw,orthometric_height)
Convert the gravity anomaly data points into a grid
Ex Code -
[Latitude_grid,Longitude_grid] = create_grid(latitude,longitude,0.01);
g_atm_and_s_mw = griddata(latitude,longitude,g_atm_and_s_mw,Latitude_grid,Longitude_grid);
-
Download the DEM data from here in ASCII format and of 5'x5' resolution. Remember to use all the tiles that cover up your intereseted area.
-
create a string array of filePaths
files = ["file1path","file2path"...]
Apply gravimetric terrain reduction to compute the Faye anomaly .
Code -
gfaye = calc_gfaye(g_atm_and_s_mw,files,Latitude_grid,Longitude_grid,0.01);
Calculate disturbing potential by using Stokes integral.
[`Tr = calc_disturb_potential(Latitude_grid,Longitude_grid,gfaye);`](https://github.com/Jainam-IITK/CE678_PROJECT/blob/main/calc_disturb_potential.m)By using Bruns’s formula, calculate undulation.
[`Nr = calc_undulation(Tr,Latitude_grid,Ellipsoid);`](https://github.com/Jainam-IITK/CE678_PROJECT/blob/main/calc_undulation.m)Restore the undulation corresponding to the removed long-wavelength gravity anomaly
of the GGM model . You will get a co-geoid
N_co_geoid = calc_n_cogeoid(Nr,nggm,latitude,longitude,Latitude_grid,Longitude_grid);
Question | PROGRESS |
---|---|
1 | Done |
2 | Done |
3 | Done |
4 | Done |
5 | Done |
6 | Done |
7 | Done |
8 | Done |
9 | Done |
10 | Done |
11 | Partially Done |