% Starting point here is that you have dataset.txt in your directory.
% the following script reads in the text file and converts it into MATLAB
% variables lat, lon, elev, freeAnomaly, bourguerAnomaly, and more.

%ProcessGravFile
load CascadiaGravity

% You should make a plot of where the gravity data lie by adding the lat's
% and lon's to a base map of the Cascadia area using the mapp function
%%
figure(99)
 args.file='puget_big';
 args.map_limits=[ -125 -121 46 49];
% args.map_limits=[ -124 -121 48 50];
 mapp(args);
 hold on
 plot(data(:,1),data(:,2),'o')
set(gca,'dataaspectratio',[1,cosd(mean(args.map_limits(3:4))),1])
 hold off
%%
% One thing to note at this point is where the gravity data are and are
% not.  You should see that although there are lots of measurement
% locations, a few spots have been missed and might make a good independent
% project for an ambitious undergrad.

% The problem now is to make a "pretty" and easy to interpret visual of the
% gravity data.  The data can be plotted as if they were a topographic "surface" -
% higher gravity sticks up and lower gravity forms the low points.

% Most computer applications can plot surfaces which are points defined on
% a regular (equally spaced) grid.  Your data are not equally spaced. 
% MATLAB has a function that interpolates between irregularly spaced
% points and return a matrix of points on a regular grid.  To do this,
% start by defining the grid:

lat_eq=39:.01:52;  % lats and lons in 0.01 degree steps
lon_eq=-115:-.01:-126;
%%
%args.map_limits=[ -123.5 -121.5 46.8 48.6];
lat_eq=args.map_limits(3):.005:args.map_limits(4);  % lats and lons in 0.01 degree steps
lon_eq=args.map_limits(2):-.005:args.map_limits(1);
grav=griddata(data(:,1),data(:,2),data(:,4),lon_eq(:)',lat_eq(:));
%%

% the MATLAB function is griddata.  The (:) put after longitudes and
% latitudes is a trick to make those variables be vertical columns. The '
% after the longitude variables is the trick to transpose longitudes into
% row vectors.  We want the maps to show longitudes on the horizontal and
% latitudes on the vertical and these steps insure that is the result.

%grav=griddata(lon(:)',lat(:),bouguerAnomaly(:),lon_eq(:)',lat_eq(:));
%%

%%
% grav is a matrix of Bouguer gravity anomaly values on an equally spaced
% grid of latitudes and longitudes.  You can now plot them with the
% following command:

mesh(lon_eq,lat_eq,grav)
%%
%an alternate function:
figure(2)
 mapp(args);
 hold on
 %plot(data(:,1),data(:,2),'o','markerfacecolor','k','markeredgecolor','k','markersize',2)
 mesh(lon_eq,lat_eq,grav-60)
%set(gca,'dataaspectratio',[1,cosd(mean(args.map_limits(3:4))),1])
 hold off



% makes the surface look different.  Both have features like "shading" that can be set to
% values other than the default.  

% your viewing perspective can be changed with the mouse (by clicking on
% the appropriate tool on the plot toolbar), or by using the command "view".
% For example looking from the south at an elevation angle of 72 degrees:

view(0,90)

% you can plot a base map over the gravity surface by "holding" the plot
% and adding the map.  However, the gravity anomaly goes above zero and obscures
% the base map so I find it necessary to shift the gravity data down by a
% constant value.
 
%%
%surf(lon_eq,lat_eq,grav-60)
figure(3)
mesh(lon_eq,lat_eq,grav-60)
 args.file='puget_big';
 args.map_limits=[ -123 -121 48 51];
 hold on
 mapp(args);
%axis equal
axis vis3d
 hold off
%%

% play around with these plots. Change your view point.  Look at them from the
% top and the side.
% See what the trends are. What are the obvious features. How do they 
% compare with geologic maps?


