Naval Postgraduate School

OC 3902

FUNDAMENTALS OF MAPPING, CHARTING, and GEODESY

LABORATORY 1

Geoid and Geoid Undulations

Objectives:

1. To compare the geoid undulation among three gravity models.

2. To understand the error sources in the WGS84-U (unclassified) used in many GPS receivers.

3. To interpolate the global 1° geoid undulation data to any interesting locations.

---

---

A. Preliminaries

1. Copy the directory lab1 with all the data files to your working director from \\lrcapps\common$\OC3902\FY09Qtr1, and work in your own directories.
2. If MATLAB is not accessible in your desktop, create a network place to \\Lrcapps\Matlab71, and access MATLAB from this network place.
3. Use the three global geoid undulation datasets which are already saved in a matlab format data file (Lab1.mat) and a global coast line data globbd.dat.:

The file Lab1.mat has five variables:

lon: Longitude

lat: Latitude

ht84: Geoid undulation height of WGS 84 model with the 180^th order (unclassified)

htosu: Geoid undulation height of OSU 360th order model    

htggm02:    Geoid undulation height of GGM02 (GRACE Gravity Model V.2) model with EMG96 

4. Interesting locations:

     ==================================== 
          Site              Location  
     ==================================== 
 
 Monterey Bay       36.5°N, -122.0° 
 
 Spit, Croatia  43.5°N,   16.5°
 
 North Korea  39.7°N,  124.5°
 
 Yellowstone Park   45.0°N, -111.0°
 
     ==================================== 
 

---

B. Procedures

---

1. Open MATLAB (Double click MATLAB 6.5). Change current directory to your oc3902/lab1 directory.

2. Load wgs84 osu Geoid and GGM02 geoid undulation data:

>> load Lab1.mat;
  
 Load global coast line data as:
 
>> load globbd.dat;
 
you can use command whos to see the detail information of the data you already have.
 >> whos

3. View 3-D Geoid undulation height using surf plot and compare the results of WGS84 OSU and GGM02 models:

Make one figure for WGS84 undulation height
 
>>  figure; surf(lon,lat,ht84);
>> shading flat;
>> hold on; drawmap3d(globbd,-150);
>> title('WGS84'); set(gca,'xlim',[0,360],'ylim',[-90,90]);
>> xlabel('Longitude'); ylabel('Latitude'); zlabel('Undulation Height(m)');
 
  make another figure for OSU undulation height
 
>> figure; surf(lon,lat,htosu);
>> shading flat;
>> hold on; drawmap3d(globbd,-150);
>> title('OSU'); set(gca,'xlim',[0,360],'ylim',[-90,90]);
>> xlabel('Longitude'); ylabel('Latitude'); zlabel('Undulation Height(m)');
 
   make another figure for GGM02 undulation height
 
>> figure; surf(lon,lat,htggm02);
>> shading flat;
>> hold on; drawmap3d(globbd,-150);
>> title('GGM02'); set(gca,'xlim',[0,360],'ylim',[-90,90]);
>> xlabel('Longitude'); ylabel('Latitude'); zlabel('Undulation Height(m)');
 
 
It is very hard to see the difference of these three data> Check the maximum, minimum, and root mean square of these undulation heights.
 
>> sprintf('Maximum: %10.2f %10.2f %10.2f\n',max(ht84(:)), max(htosu(:)), max(htggm02(:)))
>> sprintf('Minimum: %10.2f %10.2f %10.2f\n',min(ht84(:)), min(htosu(:)), min(htggm02(:)))
>> sprintf('RMS: %10.2f %10.2f %10.2f\n', rms(ht84), rms(htosu), rms(htggm02))
 
  When  the three heights are too close to see the difference from 3-D surf plot, the line plot will be better to compare the difference of the cross section. 

4. Compare the three undulation heights along some cross sections.

Latitude cross-section (Longitude: 81E)
 
>> xlon=81;
>> while(xlon<0), xlon=xlon+360; end
>> while(xlon>360), xlon=xlon-360; end
>> [mmm,ii]=min(abs(lon-xlon)); ii=ii(1);
>> figure;
>> plot(lat,ht84(:,ii)); hold on;
>> plot(lat,htosu(:,ii), 'r');
>> plot(lat,htggm02(:,ii), 'g');
>> legend('WGS84', 'OSU', 'GGM02');
>> xlabel('Latitude'); ylabel('Undulation height (m)');
 
  Longitude cross-section (Latitude: 41N)
 
>> ylat=41;
>> [mmm,jj]=min(abs(lat-ylat)); jj=jj(1);
>> figure;
>> plot(lon,ht84(jj,:)'); hold on;
>> plot(lon,htosu(jj,:)', 'r');
>> plot(lon,htggm02(jj,:)', 'g');
>> legend('WGS84', 'OSU', 'GGM02');
>> xlabel('Longitude'); ylabel('Undulation height (m)');
 

5. Perform the cross-section comparison. (wgs84htsurf.m and compwgsosu.m)

   >> close all;
   >> wgs84htsurf;
  >> complonlat;

6. Interpolate to some interesting locations:

   >> xloc=[-122,16.5,124.5,-111]; yloc=[36.5,43.5,39.7,45];
 
  Note: as 0<=lon<=360, so for interpolation, xloc must be in [0,360] region.
 
 You can convert original xloc as follow:
 
   >> ii = find(xloc<0); % find any lon <0 index
   >> xloc(ii)=xloc(ii)+360; % add 360 to lon < 0 locations only.
   >> disp(xloc); % display xloc to see if there is any <0 value.
   >> ht1=interp2(lon,lat,ht84,xloc,yloc);
   >> ht2=interp2(lon,lat,htosu,xloc,yloc);
   >> ht3=interp2(lon,lat,htggm02,xloc,yloc);
   >> disp([ht1;ht2;ht3]); 
 
  Monterey CA: Undulation height -34.9 m. 

7. Compare the three geoid undulations on a fly object trail. (trail.m)

>> trail;

Save as a jpeg file as:

>> print –djpeg lab1.jpeg;

---

Lab Results Turned-In

In your interesting area, select a fly path, compute the geoid undulations, and save the results as jpeg files (lab1a.jpeg and lab1b.jpeg). Email the two jpeg files to: fan@nps.edu

---

  (Back to "Labs List" page.)

Last Updated Sep 29, 2008
POC: Peter Chu