% PLOT_GRAV - plots reduced gravity data over Mt Soledad using griddata routine
% Steven Constable 

% assumes latitudes, longitudes, and Bouguer gravity are stored in
% longitude (degrees), latitude (degrees), and gravity (milligals), respectively

%-----------------------------------------------------------------------
% all the data files we need

topoFile = 'ljtopo1.mat';        % MAT file with topography grid 
coastFile = 'coastline.latlong';  % Latitude/Longitude of coast line
  
faults1File = 'ljfaults1.latlong'; % 2 estimates of Rose Canyon fault location
faults2File = 'ljfaults2.latlong'; % 2 column data latitude & longitude

roadsFile = 'ljroads.latlong';  % 2 columns latitude, longitude
  
benchmarkFile = 'ljbenchmarks.latlong'; % 4 columsn ID, lat, long, height

gravFile = 'processedgravity2025.txt';  % processed gravity:
% processed (latitude, free air, Bouguer) gravity data stored as
% longitude (degrees), latitude (degrees), and gravity (milligals)

radc = pi/180;       % degrees in one radian

%-----------------------------------------------------------------------
% Create a plot of gravity

grav = load(gravFile);

% create vectors of x and y coordinates using latitudes and longitudes from topo map dimensions:
x_min = -117.29;
x_max = -117.19;
x_n = 85;
xmap = linspace(x_min, x_max, x_n);

y_min = 32.78;
y_max = 32.88;
y_n = 112;
ymap = linspace(y_min, y_max, y_n);

% meshgrid creates two arrays, containing x coords and y coords for entire map
[XI,YI] = meshgrid(xmap,ymap);   

% griddata projects the unevenly spaced real data onto the entire grid
map_grid = griddata(grav(:,1),grav(:,2),grav(:,3),XI,YI); 

% now we have evenly spaced data to plot using pcolor
figure
pcolor(xmap,ymap,map_grid);
shading flat;
axis equal
axis tight

hold on
colorbar
title('Rose Canyon Fault Gravity, mgal') 
xlabel('Longitude (dec deg)')
ylabel('Latitude (dec deg)');

axis ([x_min x_max y_min y_max])  
% Ensure that the aspect ratio of the plot has equal distance for lat and long
ymid = (y_min + y_max)/2.;    % Average latitude of grid
set(gca,'DataAspectRatio',[1/cos(ymid*radc),1,1]);


%-----------------------------------------------------------------------
% add the station locations

plot(grav(:,1),grav(:,2),'.k', 'markersize',20);
% ind = [1:length(grav(:,1))]';
% text(grav(:,1),grav(:,2),num2str(ind));

                     
%-----------------------------------------------------------------------
% Plot the coast line

coastline = load(coastFile);
coast.lon = coastline(:,1);
coast.lat = coastline(:,2);

lh3 = plot(coast.lon,coast.lat,'g');
set(lh3,'LineWidth',2,'Color',[0,0.5,0]); % Make the line thicker and darker green

%-----------------------------------------------------------------------
% Add the Rose Canyon fault system

faults = load(faults1File);
ljfault.lon = faults(:,1);
ljfault.lat = faults(:,2);

faults = load(faults2File);
ljfault2.lon = faults(:,1);
ljfault2.lat = faults(:,2);

lh  = plot(ljfault.lon,ljfault.lat,'y');
set(lh,'LineWidth',2);

lh2 = plot(ljfault2.lon,ljfault2.lat,'g');
set(lh2,'LineWidth',2);

%-----------------------------------------------------------------------
% Add the city roads

ljroads = load(roadsFile);
ljroad.lon = ljroads(:,1);
ljroad.lat = ljroads(:,2);

plot(ljroad.lon,ljroad.lat,'r');



