Naval Postgraduate School OC3902
FUNDAMENTALS OF MAPPING,
CHARTING, AND GEODESY
Lab9 (GPS Lab B): GPS RESEIVER (Data
Analysis)
In this lab we will closely examine two types of GPS data, Civilian uncorrected and Differential corrected, collected over a long period of time (one day with each second one data). We will compare the 2 types of GPS data and examine their statistics. The 2 types of GPS data will be collected simultaneously at a known reference location on top of Spanagel Hall (SP2). The lab is designed to help understanding of the limitations and temporal variations of GPS.
There are three matlab format data files need to be loaded.
GPS.mat – non corrected GPS data with variables: GPS_Lon, GPS_Lat, GPS_Elv, GPS_SV, GPS_PDOP.
DGPS.mat – Differential GPS data with variables: DGPS_Lon, DGPS_Lat, DGPS_Elv, DGPS_SV, DGPS_PDOP.
SP2.mat - Spanagel Hall marker #2 reference location data with variables: SP2_Lon, SP2_Lat, SP2_Elv.
In this lab, we will plot the raw data, Quantify the Errors, and plot the distributions of these data value.
This can be done by using the ‘load’ command or by using the Open function in the menu pull down. If using the load command format is as follows: load ‘Your selected directory name\DGPS’
: load ‘Your selected directory name\GPS’
: load ‘Your selected directory name\SP2’
Or when using the Open function navigate to the files in question and select. This must be done for each file.
The variables that are now loaded (you can see them either in your workspace window, or by using the command ‘whos’ are of the following format. Part1_Part2
The Part1 before the underscore tells you what data source it is from. DGPS = Differential GPS ,
GPS = non corrected GPS,
SP2= Spanagel Hall marker #2 reference location.
The Part2 after the underscore tells what type of data it is.
Lat = Latitude in decimal Degrees
Lon= Longitude in decimal Degrees
Elv = Elevation in Meters
SV = number of Satellite vehicles used in the position calculation
PDOP= Precision Dilution of Position
You will note that DGPS and GPS data has multiple data points. Each data point was taken every second.
1.
Plot Raw
Data (rawdata.m)
In order to compare the GPS and DGPS results plot both GPS (red)and DGPS (blue) raw data together. There are three major variables Longitude, Latitude and Elevation will be plot in one paper. Let the top plot will be Longitude/time, the middle one will be Latitude/time and the bottom one will be Elevation/time. The units of time is better to be hour and as they have same time, as only bottom one need xlabel. The follow program (GPS_plot.m) plot the top one (Longitude/time) raw data.
time=[1:length(GPS_Lon)]/3600; % convert time from second to hour
Glonmean=mean(GPS_Lon); Glatmean=mean(GPS_Lat);
Gelvmean=mean(GPS_Elv);
DGlonmean=mean(DGPS_Lon); DGlatmean=mean(DGPS_Lat);
DGelvmean=mean(DGPS_Elv);
yl=0.08; dy=0.28; ddy=dy+0.02;
figure('units','inches','position',[1,2,5,7]);
axes('position',[0.15,yl+2*ddy,0.82,dy],'xlim',[0,24],'xtick',[4:4:20],'box','on','Fontsize',8);
hold on;
set(gca,'xticklabel',[]);
plot(time,GPS_Lon,'r');
plot([0,24],Glonmean*[1,1],'r');
plot(time,DGPS_Lon,'b');
plot([0,24],DGlonmean*[1,1],'b');
ylm=get(gca,'ylim');
dym=0.1*diff(ylm);
y1=ylm(2)-dym; y2=ylm(1)+dym;
text(12,y1,['Mean GPS Lon: ',num2str(Glonmean,8)],'HorizontalAlignment','center','Fontsize',8);
text(12,y2,['Mean DGPS Lon: ',num2str(DGlonmean,8)],'HorizontalAlignment','center','Fontsize',8);
ylabel('Longitude, (Degrees)');
title('Red is GPS, Blue is DGPS')
Please finish the Latitude/time and Elevation/time raw data plots as middle and bottom subplot. On the bottom one make an xlabel. Save the figure as a jpg file.
2.
Quantify the
Errors and correlation of Number of Satellite Vehicles (errplot.m)
a. Convert to Meters:
While the previous plot is nice looking, it would more desireable if the positions were displayed as a difference from a reference and in meters instead of geodetic degrees. Recall the basic conversions for geodetic degrees to meters is a constant for latitude but varies for longitude.
Conversion numbers needed to get your conversion factor.
1 degree = 60 nautical miles
1 nautical mile = 1852 meters
Make a variable called ‘conv_factor’(=60*1852) and assign a value that will allow you to multiply it by your GPS latitude variables and get meters.
erGPS_LatM= conv_factor*(GPS_Lat – SP2_Lat);
erGPS_LonM=conv_factor.* cos((pi/180).*GPS_Lat).*(GPS_Lon– SP2_Lon);
b. Make Plots
Lets plot correlation of # of satellites to error as the three rows and two columns plot. The first column is Root Mean Square Error, the second column is Standard Deviation, and tht three rows will respond to Longitude (X), Latitude (Y) and Elevation (Z). The Total 3-D RMSE and STD will be written in the title. The follow program (compsv.m) plots the first row (X) and the all calculation.
fact=60*1852;
factx=fact*cos(SP2_Lat*pi/180);
erGPS_LonM=(GPS_Lon-SP2_Lon)*factx;
erDGPS_LonM=(DGPS_Lon-SP2_Lon)*factx;
erGPS_LatM=(GPS_Lat-SP2_Lat)*fact;
erDGPS_LatM=(DGPS_Lat-SP2_Lat)*fact;
erGPS_Elv=GPS_Elv-SP2_Elv;
erDGPS_Elv=DGPS_Elv-SP2_Elv;
gsvmin=min(GPS_SV);
gsvmax=max(GPS_SV);
dgsvmin=min(DGPS_SV);
dgsvmax=max(DGPS_SV);
GPS3der=num2str(sqrt(mean(erGPS_LonM.^2+erGPS_LatM.^2+erGPS_Elv.^2)),3);
DGPS3der=num2str(sqrt(mean(erDGPS_LonM.^2+erDGPS_LatM.^2+erDGPS_Elv.^2)),3);
GPS3dstd=num2str(std(sqrt(erGPS_LonM.^2+erGPS_LatM.^2+erGPS_Elv.^2)),3);
DGPS3dstd=num2str(std(sqrt(erDGPS_LonM.^2+erDGPS_LatM.^2+erDGPS_Elv.^2)),3);
for
i=4:gsvmax
[j]=find(GPS_SV==i);
if ~isempty(j)
erlon(i-3)=sqrt(mean(erGPS_LonM(j).^2));
erlat(i-3)=sqrt(mean(erGPS_LatM(j).^2));
erelv(i-3)=sqrt(mean(erGPS_Elv(j).^2));
stdlon(i-3)=std(GPS_Lon(j))*factx;
stdlat(i-3)=std(GPS_Lat(j))*fact;
stdelv(i-3)=std(GPS_Elv(j));
end
end
for
i=4:dgsvmax
[j]=find(DGPS_SV==i);
if ~isempty(j)
derlon(i-3)=sqrt(mean(erDGPS_LonM(j).^2));
derlat(i-3)=sqrt(mean(erDGPS_LatM(j).^2));
derelv(i-3)=sqrt(mean(erDGPS_Elv(j).^2));
dstdlon(i-3)=std(DGPS_Lon(j))*factx;
dstdlat(i-3)=std(DGPS_Lat(j))*fact;
dstdelv(i-3)=std(DGPS_Elv(j));
end
end
n1=[4:gsvmax];
n2=[4:dgsvmax];
figure('units','inches','position',[4,2,5,7]);
xl=0.1;
dx=0.41; ddx=dx+0.065; yl=0.06; dy=0.27; ddy=dy+0.025;
axes('position',[xl,yl+2*ddy,dx,dy],'xlim',[4,12],'box','on','fontsize',8);
hold
on; set(gca,'xticklabel',[]);
title(str2mat('RMSE',['GPS(red) Total3DRMSE(m):',GPS3der],
['DGPS(blue) Total3DRMSE(m):',DGPS3der]));
plot(n1,erlon,'r*-');
plot(n2,derlon,'b*-');
ylabel('X (m)');
axes('position',[xl+ddx,yl+2*ddy,dx,dy],'xlim',[4,12],'box','on','fontsize',8);
hold
on; set(gca,'xticklabel',[]);
title('STD
(Red is GPS, Blue is DGPS)');
title(str2mat('STD',['GPS(red) Total3DSTD(m):',GPS3dstd],
['DGPS(blue) Total3DSTD(m):',DGPS3dstd]));
plot(n1,stdlon,'r*-');
plot(n2,dstdlon,'b*-');
Please finish the rest rows (Y for middle row and Z for bottom row). Make xlabel on the bottom row. Save as jpg file.
3.
Distribution
Plot (Histogram – histplot.m)
Data distribution is also important in data analysis. This plot will also be two columns and three rows. The first column is GPS, and the second column is DGPS. The three rows respond to Longitude, Latitude and Elevation. The first row plot can be as the follow (comphist.m).
figure('units','inches','position',[7,2,5,7]);
xl=0.12;
dx=0.41; ddx=dx+0.045; yl=0.06; dy=0.25; ddy=dy+0.07;
axes('position',[xl,yl+2*ddy,dx,dy],'xtick',[],'box','on','fontsize',8);
hold
on;
title('Distribution
of GPS');
[N,X]=hist(GPS_Lon,20);
hist(GPS_Lon,20);
ylm=get(gca,'ylim');
plot(SP2_Lon*[1,1],ylm,'r');
text(SP2_Lon,ylm(2),num2str(SP2_Lon,10),'color','r','fontsize',8,...
'horizontalAlignment','right','verticalalignment','bottom','rotation',90);
ii=find(N>0.02*max(N));
X=X(ii([1,end]));
plot([1;1]*X,ylm'*[1,1],'--');
ym=mean(ylm);
text(X(1),ym,num2str(X(1),10),'fontsize',8,'rotation',90,...
'horizontalAlignment','center','verticalalignment','bottom');
text(X(2),ym,num2str(X(2),10),'fontsize',8,'rotation',90,...
'horizontalAlignment','center','verticalalignment','top');
xlabel('Longitude
(degrees)'); ylabel('Frequency');
axes('position',[xl+ddx,yl+2*ddy,dx,dy],'xtick',[],'box','on','fontsize',8);
hold
on;
title('Distribution
of DGPS');
[N,X]=hist(DGPS_Lon,20);
hist(DGPS_Lon,20);
ylm=get(gca,'ylim');
plot(SP2_Lon*[1,1],ylm,'r');
text(SP2_Lon,ylm(2),num2str(SP2_Lon,10),'color','r','fontsize',8,...
'horizontalAlignment','right','verticalalignment','bottom','rotation',90);
ii=find(N>0.02*max(N));
X=X(ii([1,end]));
plot([1;1]*X,ylm'*[1,1],'--');
ym=mean(ylm);
text(X(1),ym,num2str(X(1),10),'fontsize',8,'rotation',90,...
'horizontalAlignment','center','verticalalignment','bottom');
text(X(2),ym,num2str(X(2),10),'fontsize',8,'rotation',90,...
'horizontalAlignment','center','verticalalignment','top');
xlabel('Longitude
(degrees)');
Please finish the middle and bottom rows (Latitude and Elevation). Note: as these plots have different x variable, so the xlabel should be made individual. Save as a jpg file.
Lab Results Turned in:
Email these three jpg files to fan@nps.edu
(Back to "Labs
List" page.)
Last Updated 29 Sep. 2008
POC: Peter Chu