% PLOT_MTSOLEDAD_MAP - plots a map of La Jolla/Mt Soledad area topography adding
% roads, Rose canyon fault, city benchmarks and the locations of the class gravity measurements

% Alistiar Harding, with later mods by Steven Constable.

%-----------------------------------------------------------------------
% 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

gravdataFile = 'gravity2025.txt';  % class gravity database:
% 11 cols: No, ID, Lat, Long, height, northing, easting, hours, mins, meter reading, raw gravity
% gravity is corrected for meter calibration and drift only
% heights are in feet
% northing and easting are in 100's of feet

radc = pi/180;       % degrees in one radian

%-----------------------------------------------------------------------
% Create a plot of topography

load(topoFile)       % Load ljh - matrix with la jolla topography
                     %  Structures defining grid coords
                     % ljgeo   - coordinate parameters in (global) lat/lon
                     % ljlocal - coordinate parameter  in local km
                     
ljlon = linspace(ljgeo.x_min,ljgeo.x_max,ljgeo.nx);   % create coordinate vectors
ljlat = linspace(ljgeo.y_min,ljgeo.y_max,ljgeo.ny);

% Create a filled contour plot of topography

figure
contour_levels = [1, 10:10:240];         
contourf(ljlon, ljlat,ljh,contour_levels); % Create a filled contour map
hold on
% axis ([-117.29 -117.2 32.78 32.88])  
axis ([-117.29 -117.18 32.78 32.90])  

colorbar
title('Rose Canyon Fault + topography, m') 
xlabel('Longitude (dec deg)')
ylabel('Latitude (dec deg)');

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


% Install a customized colormap
mycmap = copper;    % Matlab has some standard colormaps, help copper gives list 
colormap(mycmap);
brighten(0.7);        % Make the colors in the colormap brighter
mycmap = colormap;
mycmap(1,:) = mycmap(2,:)*0.9;  % Replace the too dark 1st color entry
colormap(mycmap);

%-----------------------------------------------------------------------
% 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');


%-----------------------------------------------------------------------
% Add positions of the city benchmarks

ljbenchmarks = load(benchmarkFile);

ft2m = 1/3.2808399;   % convert ft to meters

ljbench.id  = ljbenchmarks(:,1);
ljbench.lon = ljbenchmarks(:,2);    % dec deg
ljbench.lat = ljbenchmarks(:,3);    % dec deg
ljbench.ht  = ljbenchmarks(:,4) *ft2m;     % ht in meters
plot(ljbench.lon,ljbench.lat,'c.');

%-----------------------------------------------------------------------
% Add locations where the class collected data

gravdata = load (gravdataFile);
gravlon = gravdata(:,4);
gravlat = gravdata(:,3);

plot(gravlon,gravlat,'ow','markerfacecolor','w','markersize',6);

plot(gravlon(225:231),gravlat(225:231),'om','markersize',7);  % 2025 data

plot(gravlon(35:48),gravlat(35:48),'ok','markersize',7);  % Nettleton's data

legend('topo', 'coastline', 'Faults1', 'Faults2','roads', 'benchmarks', 'gravity sites', '2025 data', 'Nettleton data','Location','NorthWest')

axis([-117.285 -117.19 32.78 32.88])
