From f4e3169bd19dcfd246ac4bc8d06e72ace6eacadc Mon Sep 17 00:00:00 2001
From: Chris Leaman <ckleaman@gmail.com>
Date: Tue, 20 Nov 2018 10:48:27 +1100
Subject: [PATCH] Update beach orientations MATALB to run on command

---
 src/data/beach_orientations.m | 76 +++++++++++++++++++++++++++++------
 1 file changed, 64 insertions(+), 12 deletions(-)

diff --git a/src/data/beach_orientations.m b/src/data/beach_orientations.m
index d537ad6..df0a9a5 100644
--- a/src/data/beach_orientations.m
+++ b/src/data/beach_orientations.m
@@ -1,19 +1,39 @@
-% Calculate orientation the beach profile at each unique site and save to .mat file
+% Calculate orientation the beach profile at each unique site and save to .mat file. Orientation is
+% the number of degrees, anticlockwise from east, perpendicular to the shoreline (pointing towards
+% land).
+
+%% Setup
 % Needs the following coastal tools:
-% J:\Coastal\Tools\MALT Logspiral Transformation
-% J:\Coastal\Tools\Coordinate Transformations
+addpath('J:\Coastal\Tools\MALT Logspiral Transformation')
+addpath('J:\Coastal\Tools\Coordinate Transformations')
 
 clear 
 clc
 
+%% Options
+% Where is the profiles file located? This should contain a structure including the .lat and .lon
+% for each analysed cross section
+profilesFile = '..\..\data\raw\processed_shorelines\profiles.mat';
+
+% Where should we store the processed beach orientations?
+outputFile = '..\..\data\raw\processed_shorelines\orientations.mat';
+
+% How far in meters does the profile extend towards land and sea? Used to provide end points of the
+% cross section
+distance = 200;
+
+
+%% Script
+
 % Load profile data, this is where we want to calculate orientations.
 warning('off','all')
-data = load('C:\Users\z5189959\Desktop\nsw_2016_storm_impact\data\raw\processed_shorelines\profiles.mat');
+data = load(profilesFile);
 data = data.data;
 
+% Save results to variable
 output = [];
 
-for ii = 1:n
+for ii = 1:length(data)
   disp(num2str(ii))
   lat = data(ii).lat;
   lon = data(ii).lon;
@@ -21,14 +41,41 @@ for ii = 1:n
   [x,y,utmzone] = deg2utm(lat,lon);
 
   if strcmp(beach, 'BOOM') == 1 || strcmp(beach, 'HARGn') == 1 || strcmp(beach, 'BILG') == 1 || strcmp(beach, 'HARGs') == 1 || strcmp(beach, 'DEEWHYn') == 1
-    % log spiral transformation file is out of date. Talk to Mitch
-    continue
-  end
+    % These are straight beaches, load the transformation file directly and read the rotation angle.
+    parameterDir = 'J:\Coastal\Tools\MALT Logspiral Transformation';
+    parameterFile = [parameterDir, filesep, beach, '.mat'];
+    parameterMat = load(parameterFile);
+    fields = fieldnames(parameterMat);
+    field = fields(1,1);
+    site = getfield(parameterMat, field{1});
+    rot_angle = site.rot_angle;  % Angle of the shoreline counter clockwise from east
+    
+    % Figure out end points in utm coordinates
+
+    x_land = x - distance * cos(deg2rad(rot_angle));
+    y_land = y + distance * sin(deg2rad(rot_angle));
+    x_sea = x + distance * cos(deg2rad(rot_angle));
+    y_sea = y - distance * sin(deg2rad(rot_angle));
+    
+    [lat_land,lon_land] = utm2deg(x_land,y_land,utmzone);
+    [lat_sea,lon_sea] = utm2deg(x_land,y_land,utmzone);
     
-  if strcmp(beach, 'AVOCAs') == 1
-    % negative solution. Talk to Mitch
+    row.lat_center = lat;
+    row.lon_center = lon;
+    row.lat_land = lat_land;
+    row.lon_land = lat_land;
+    row.lat_sea = lat_sea;
+    row.lon_sea = lon_sea;
+    row.orientation = rot_angle + 90;  % Tangent to shoreline towards line 
+    row.beach = beach;
+    output = [output; row];
     continue
   end
+    
+%   if strcmp(beach, 'AVOCAs') == 1
+%     % negative solution. Talk to Mitch
+%     continue
+%   end
   
   % Get the sp log spiral transformed coordinates
   xyz.x = x;
@@ -56,7 +103,12 @@ for ii = 1:n
   [lat_sea,lon_sea] = utm2deg(xyz_sea.x,xyz_sea.y,utmzone);
 
   % Orientation in degrees anticlockwise from east, pointing towards land
-  orientation = radtodeg(atan2((xyz_land.y - xyz_sea.y), (xyz_land.x - xyz_sea.x)));
+  try
+    orientation = radtodeg(atan2((xyz_land.y - xyz_sea.y), (xyz_land.x - xyz_sea.x)));
+  catch
+    disp(['Cannot calculate orientation: ' beach])
+    continue
+  end
   
   row.lat_center = lat;
   row.lon_center = lon;
@@ -70,4 +122,4 @@ for ii = 1:n
 end
 warning('on','all')
 
-save('orientations.mat','output','-v7')
+save(outputFile','output','-v7')