CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
create_profiles.cpp
Go to the documentation of this file.
1
12
13/* ==============================================================================================================================
14
15 This file is part of CoastalME, the Coastal Modelling Environment.
16
17 CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
20
21 You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22
23==============================================================================================================================*/
24#include <assert.h>
25
26#include <cstdio>
27#include <cmath>
28#include <cfloat>
29
30#include <iostream>
31using std::cerr;
32// using std::cout;
33using std::endl;
34using std::ios;
35
36#include <algorithm>
37using std::find;
38using std::sort;
39
40#include <utility>
41using std::make_pair;
42using std::pair;
43
44#include <random>
45using std::normal_distribution;
46
47#include "cme.h"
48#include "simulation.h"
49#include "coast.h"
50#include "2d_point.h"
51#include "2di_point.h"
52
53namespace
54{
55//===============================================================================================================================
57//===============================================================================================================================
58bool bCurvaturePairCompareDescending(const pair<int, double>& prLeft, const pair<int, double>& prRight)
59{
60 // Sort in descending order (i.e. most concave first)
61 return prLeft.second > prRight.second;
62}
63} // namespace
64
65//===============================================================================================================================
67//===============================================================================================================================
69{
71 LogStream << m_ulIter << ": Creating profiles" << endl;
72
73 for (unsigned int nCoast = 0; nCoast < m_VCoast.size(); nCoast++)
74 {
75 int nProfile = 0;
76 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
77
78 // Create a bool vector to mark coast points which have been searched
79 vector<bool> bVCoastPointDone(nCoastSize, false);
80
81 // Now create a vector of pairs: the first value of the pair is the coastline point, the second is the coastline's curvature at that point
82 vector<pair<int, double>> prVCurvature;
83
84 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
85 {
86 double dCurvature;
87
88 if (m_VCoast[nCoast].pGetCoastLandform(nCoastPoint)->nGetLandFormCategory() != LF_CAT_INTERVENTION)
89 {
90 // Not an intervention coast point, so store the smoothed curvature
91 dCurvature = m_VCoast[nCoast].dGetSmoothCurvature(nCoastPoint);
92 }
93
94 else
95 {
96 // This is an intervention coast point, which is likely to have some sharp angles. So store the detailed curvature
97 dCurvature = m_VCoast[nCoast].dGetDetailedCurvature(nCoastPoint);
98 }
99
100 prVCurvature.push_back(make_pair(nCoastPoint, dCurvature));
101 }
102
103 // Sort this pair vector in descending order, so that the most convex curvature points are first
104 sort(prVCurvature.begin(), prVCurvature.end(), bCurvaturePairCompareDescending);
105
106 // // DEBUG CODE =======================================================================================================================
107 // for (int n = 0; n < prVCurvature.size(); n++)
108 // {
109 // LogStream << prVCurvature[n].first << "\t" << prVCurvature[n].second << endl;
110 // }
111 // LogStream << endl << endl;
112 // // DEBUG CODE =======================================================================================================================
113
114 // And mark points at and near the start and end of the coastline so that they don't get searched (will be creating 'special' start- and end-of-coast profiles at these end points later)
115 for (int n = 0; n < m_nCoastNormalSpacing; n++)
116 {
117 if (n < nCoastSize)
118 bVCoastPointDone[n] = true;
119
120 int const m = nCoastSize - n - 1;
121
122 if (m >= 0)
123 bVCoastPointDone[m] = true;
124 }
125
126 // Now locate the start points for all coastline-normal profiles (except the grid-edge ones), at points of maximum convexity. Then create the profiles
127 LocateAndCreateProfiles(nCoast, nProfile, &bVCoastPointDone, &prVCurvature);
128
129 // Did we fail to create any normal profiles? If so, quit
130 if (nProfile < 0)
131 {
132 string strErr = ERR + "timestep " + strDblToStr(m_ulIter) + ": could not create profiles for coastline " + strDblToStr(nCoast);
133
134 if (m_ulIter == 1)
135 strErr += ". Check the SWL";
136
137 strErr += "\n";
138
139 cerr << strErr;
140 LogStream << strErr;
141
143 }
144
145 // Locate and create a 'special' profile at the grid edge, first at the beginning of the coastline. Then put this onto the raster grid
146 int nRet = nLocateAndCreateGridEdgeProfile(true, nCoast, nProfile);
147
148 if (nRet != RTN_OK)
149 return nRet;
150
151 // Locate a second 'special' profile at the grid edge, this time at end of the coastline. Then put this onto the raster grid
152 nRet = nLocateAndCreateGridEdgeProfile(false, nCoast, ++nProfile);
153
154 if (nRet != RTN_OK)
155 return nRet;
156
157 // Insert pointers to profiles at coastline points in the profile-all-coastpoint index
158 m_VCoast[nCoast].InsertProfilesInProfileCoastPointIndex();
159
160 // // DEBUG CODE ===================================================================================================
161 // LogStream << endl << "===========================================================================================" << endl;
162 // LogStream << "PROFILES BEFORE ADDING BEFORE- AND AFTER-PROFILE NUMBERS" << endl;
163 // int nNumProfiles = m_VCoast[nCoast].nGetNumProfiles();
164 // for (int nn = 0; nn < nNumProfiles; nn++)
165 // {
166 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(nn);
167 //
168 // LogStream << nn << " nCoastID = " << pProfile->nGetProfileID() << " nGlobalID = " << pProfile->nGetProfileID() << " nGetCoastPoint = " << pProfile->nGetCoastPoint() << " pGetUpCoastAdjacentProfile = " << pProfile->pGetUpCoastAdjacentProfile() << " pGetDownCoastAdjacentProfile = " << pProfile->pGetDownCoastAdjacentProfile() << endl;
169 // }
170 // LogStream << "===================================================================================================" << endl << endl;
171 // // DEBUG CODE ===================================================================================================
172
173 // int nLastProfile;
174 // int nThisProfile;
175 CGeomProfile* pLastProfile;
176 CGeomProfile* pThisProfile;
177
178 // Go along the coastline and give each profile the number of the adjacent up-coast profile and the adjacent down-coast profile
179 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
180 {
181 if (m_VCoast[nCoast].bIsProfileAtCoastPoint(nCoastPoint))
182 {
183 pThisProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(nCoastPoint);
184 // nThisProfile = pThisProfile->nGetProfileID();
185
186 if (nCoastPoint == 0)
187 {
188 pThisProfile->SetUpCoastAdjacentProfile(NULL);
189
190 // LogStream << "nCoastPoint = " << nCoastPoint << " ThisProfile = " << nThisProfile << " ThisProfile UpCoast = " << pThisProfile->pGetUpCoastAdjacentProfile() << " ThisProfile DownCoast = " << pThisProfile->pGetDownCoastAdjacentProfile() << endl;
191
192 pLastProfile = pThisProfile;
193 // nLastProfile = nThisProfile;
194 continue;
195 }
196
197 // nLastProfile = pLastProfile->nGetProfileID();
198 pLastProfile->SetDownCoastAdjacentProfile(pThisProfile);
199 pThisProfile->SetUpCoastAdjacentProfile(pLastProfile);
200
201 if (nCoastPoint == nCoastSize - 1)
202 pThisProfile->SetDownCoastAdjacentProfile(NULL);
203
204 // LogStream << "nCoastPoint = " << nCoastPoint << " LastProfile = " << nLastProfile << " LastProfile UpCoast = " << pLastProfile->pGetUpCoastAdjacentProfile() << " LastProfile DownCoast = " << pLastProfile->pGetDownCoastAdjacentProfile() << " ThisProfile = " << nThisProfile << " ThisProfile UpCoast = " << pThisProfile->pGetUpCoastAdjacentProfile() << " ThisProfile DownCoast = " << pThisProfile->pGetDownCoastAdjacentProfile() << endl;
205
206 pLastProfile = pThisProfile;
207 // nLastProfile = nThisProfile;
208 }
209 }
210
211 // And create an index to this coast's profiles in along-coastline sequence
212 m_VCoast[nCoast].CreateProfileDownCoastIndex();
213
214 // // DEBUG CODE =======================================================================================================================
215 // for (int n = 0; n < m_VCoast[nCoast].nGetNumProfiles(); n++)
216 // {
217 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(n);
218 // CGeomProfile* pUpCoastProfile = pProfile->pGetUpCoastAdjacentProfile();
219 // CGeomProfile* pDownCoastProfile = pProfile->pGetDownCoastAdjacentProfile();
220 // int nUpCoastProfile = INT_NODATA;
221 // int nDownCoastProfile = INT_NODATA;
222 // if (pUpCoastProfile != 0)
223 // nUpCoastProfile = pUpCoastProfile->nGetProfileID();
224 // if (pDownCoastProfile != 0)
225 // nDownCoastProfile = pDownCoastProfile->nGetProfileID();
226 // LogStream << "nCoastID = " << pProfile->nGetProfileID() << "\t up-coast profile = " << nUpCoastProfile << "\t down-coast profile = " << nDownCoastProfile << endl;
227 // }
228 // LogStream << endl;
229 // // DEBUG CODE =======================================================================================================================
230
231 // // DEBUG CODE =======================================================================================================================
232 // int nProf = 0;
233 // for (int n = 0; n < nCoastSize; n++)
234 // {
235 // // LogStream << n << "\t";
236 //
237 // // LogStream << m_VCoast[nCoast].dGetDetailedCurvature(n) << "\t";
238 // //
239 // // LogStream << m_VCoast[nCoast].dGetSmoothCurvature(n) << "\t";
240 // //
241 // // if (m_VCoast[nCoast].pGetCoastLandform(n)->nGetLandFormCategory() == LF_CAT_INTERVENTION)
242 // // LogStream << "I\t";
243 //
244 // if (m_VCoast[nCoast].bIsProfileAtCoastPoint(n))
245 // {
246 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(n);
247 //
248 // LogStream << "profile " << pProfile->nGetProfileID() << " at coast point " << n << " adjacent up-coast profile = " << pProfile->pGetUpCoastAdjacentProfile() << " adjacent down-coast profile = " << pProfile->pGetDownCoastAdjacentProfile() << endl;
249 //
250 // nProf++;
251 // }
252 // }
253 // LogStream << endl;
254 // LogStream << "nProf = " << nProf << endl;
255 // // DEBUG CODE =======================================================================================================================
256
257 // // DEBUG CODE =======================================================================================================================
258 // LogStream << "=====================" << endl;
259 // for (int n = 0; n < m_VCoast[nCoast].nGetNumProfiles(); n++)
260 // {
261 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(n);
262 // int nStartPoint = pProfile->nGetCoastPoint();
263 //
264 // LogStream << n << "\t nCoastID = " << pProfile->nGetProfileID() << "\tnStartPoint = " << nStartPoint << endl;
265 // }
266 // LogStream << endl;
267 // LogStream << "=====================" << endl;
268 // // DEBUG CODE =======================================================================================================================
269 }
270
271 return RTN_OK;
272}
273
274//===============================================================================================================================
276//===============================================================================================================================
277void CSimulation::LocateAndCreateProfiles(int const nCoast, int& nProfile, vector<bool>* pbVCoastPointDone, vector<pair<int, double>> const* prVCurvature)
278{
279 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
280
281 // Work along the vector of curvature pairs starting at the convex end
282 for (int n = nCoastSize - 1; n >= 0;
283 n--)
284 {
285 // Have we searched all the coastline points?
286 int nStillToSearch = 0;
287
288 for (int m = 0; m < nCoastSize; m++)
289 if (! pbVCoastPointDone->at(m))
290 nStillToSearch++;
291
292 if (nStillToSearch == 0)
293 // OK we are done here
294 return;
295
296 // This convex point on the coastline is a potential location for a normal
297 int const nNormalPoint = prVCurvature->at(n).first;
298
299 // Ignore each end of the coastline
300 if ((nNormalPoint == 0) || (nNormalPoint == nCoastSize - 1))
301 continue;
302
303 // TODO 089 When choosing locations for profiles, do coast first then interventions
304
305 if (! pbVCoastPointDone->at(nNormalPoint))
306 {
307 // We have not already searched this coast point. Is it an intervention coast point?
308 bool bIntervention = false;
309
310 if (m_VCoast[nCoast].pGetCoastLandform(nNormalPoint)->nGetLandFormCategory() == LF_CAT_INTERVENTION)
311 {
312 // It is an intervention
313 bIntervention = true;
314 }
315
316 CGeom2DIPoint const PtiThis = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nNormalPoint);
317
318 // Create a profile here
319 int const nRet = nCreateProfile(nCoast, nCoastSize, nNormalPoint, nProfile, bIntervention, &PtiThis);
320
321 // // DEBUG CODE =================
322 // LogStream << "After nCreateProfile() ===========" << endl;
323 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(nProfile);
324 // LogStream << pProfile->nGetProfileID() << "\t";
325 //
326 // int nPointsInProfile = pProfile->nGetProfileSize();
327 //
328 // for (int nPoint = 0; nPoint < nPointsInProfile; nPoint++)
329 // {
330 // CGeom2DPoint Pt = *pProfile->pPtGetPointInProfile(nPoint);
331 // LogStream << " {" << Pt.dGetX() << ", " << Pt.dGetY() << "}";
332 // }
333 // LogStream << endl << "===========" << endl;
334 // // DEBUG CODE =================
335
336 // Mark this coast point as searched
337 pbVCoastPointDone->at(nNormalPoint) = true;
338
339 if (nRet != RTN_OK)
340 {
341 // This potential profile is no good (has hit coast, or hit dry land, etc.) so forget about it
342 // LogStream << "Profile is no good" << endl;
343 continue;
344 }
345
346 // // DEBUG CODE ===================================================================================================
347 // LogStream << endl << "===========================================================================================" << endl;
348 // LogStream << "PROFILES JUST AFTER CREATION" << endl;
349 // int nNumProfiles = m_VCoast[nCoast].nGetNumProfiles();
350 // for (int nn = 0; nn < nNumProfiles; nn++)
351 // {
352 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(nn);
353 //
354 // LogStream << nn << " nCoastID = " << pProfile->nGetProfileID() << " nGlobalID = " << pProfile->nGetProfileID() << " nGetCoastPoint = " << pProfile->nGetCoastPoint() << " pGetUpCoastAdjacentProfile = " << pProfile->pGetUpCoastAdjacentProfile() << " pGetDownCoastAdjacentProfile = " << pProfile->pGetDownCoastAdjacentProfile() << endl;
355 // }
356 // LogStream << "===================================================================================================" << endl << endl;
357 // // DEBUG CODE ===================================================================================================
358 //
359 // // DEBUG CODE ===================================================================================================
360 // LogStream << "++++++++++++++++++++++" << endl;
361 // LogStream << endl << "Just created profile " << nProfile << endl;
362 // int nProf = 0;
363 // for (int nnn = 0; nnn < nCoastSize; nnn++)
364 // {
365 // if (m_VCoast[nCoast].bIsProfileAtCoastPoint(nnn))
366 // {
367 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(nnn);
368 //
369 // LogStream << "profile " << pProfile->nGetProfileID() << " at coast point " << nnn << " adjacent up-coast profile = " << pProfile->pGetUpCoastAdjacentProfile() << " adjacent down-coast profile = " << pProfile->pGetDownCoastAdjacentProfile() << endl;
370 //
371 // nProf++;
372 // }
373 // }
374 // LogStream << endl;
375 // LogStream << "nProf = " << nProf << endl;
376 // LogStream << "++++++++++++++++++++++" << endl;
377 // // DEBUG CODE ===================================================================================================
378
379 // CGeom2DPoint PtThis = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nNormalPoint);
380 // if (m_nLogFileDetail >= LOG_FILE_ALL)
381 // LogStream << m_ulIter << ": Coast " << nCoast << " profile " << nProfile << " created at coast point " << nNormalPoint << " [" << PtiThis.nGetX() << "][" << PtiThis.nGetY() << "] = {" << PtThis.dGetX() << ", " << PtThis.dGetY() << "} (smoothed curvature = " << m_VCoast[nCoast].dGetSmoothCurvature(nNormalPoint) << ", detailed curvature = " << m_VCoast[nCoast].dGetDetailedCurvature(nNormalPoint) << ")" << endl;
382
383 // // DEBUG CODE =================================================================================
384 // if (m_pRasterGrid->m_Cell[PtiThis.nGetX()][PtiThis.nGetY()].bIsCoastline())
385 // LogStream << m_ulIter << ": cell[" << PtiThis.nGetX() << "][" << PtiThis.nGetY() << "] IS coastline, coast number = " << m_pRasterGrid->m_Cell[PtiThis.nGetX()][PtiThis.nGetY()].nGetCoastline() << endl;
386 // else
387 // LogStream << m_ulIter << ": ******* cell[" << PtiThis.nGetX() << "][" << PtiThis.nGetY() << "] IS NOT coastline" << endl;
388 // // DEBUG CODE =================================================================================
389
390 // This profile is fine
391 nProfile++;
392
393 // We need to mark points on either side of this profile so that we don't get profiles which are too close together. However best-placed profiles on narrow intervention structures may need to be quite closes
394 double dNumToMark = m_nCoastNormalSpacing;
395
396 if (bIntervention)
398
399 // If we have a random factor for profile spacing, then modify the profile spacing
401 {
402 // Draw a sample from the unit normal distribution using random number generator 0
403 double const dRand = m_dGetFromUnitNormalDist(m_Rand[0]);
404
405 double const dTmp = dRand * m_dCoastNormalRandSpacingFactor * dNumToMark;
406 dNumToMark += dTmp;
407
408 // Make sure number to mark is not too small or too big TODO 011
409 if (bIntervention)
410 {
411 dNumToMark = tMin(dNumToMark, m_nCoastNormalInterventionSpacing * 0.75);
412 dNumToMark = tMax(dNumToMark, m_nCoastNormalInterventionSpacing * 1.25);
413 }
414
415 else
416 {
417 dNumToMark = tMin(dNumToMark, m_nCoastNormalSpacing * 0.75);
418 dNumToMark = tMax(dNumToMark, m_nCoastNormalSpacing * 1.25);
419 }
420
421 // TODO 014 Assume that the above is the profile spacing on straight bits of coast. Try gradually increasing the profile spacing with increasing concavity, and decreasing the profile spacing with increasing convexity. Could use a Michaelis-Menten S-curve relationship for this i.e.
422 // double fReN = pow(NowCell[nX][nY].dGetReynolds(m_dNu), m_dDepN);
423 // double fC1 = m_dC1Laminar - ((m_dC1Diff * fReN) / (fReN + m_dReMidN));
424 }
425
426 // Mark points on either side of the profile
427 for (int m = 1; m < dNumToMark; m++)
428 {
429 int nTmpPoint = nNormalPoint + m;
430
431 if (nTmpPoint < nCoastSize)
432 pbVCoastPointDone->at(nTmpPoint) = true;
433
434 nTmpPoint = nNormalPoint - m;
435
436 if (nTmpPoint >= 0)
437 pbVCoastPointDone->at(nTmpPoint) = true;
438 }
439 }
440 }
441}
442
443//===============================================================================================================================
445//===============================================================================================================================
446int CSimulation::nCreateProfile(int const nCoast, int const nCoastSize, int const nProfileStartPoint, int const nProfile, bool const bIntervention, CGeom2DIPoint const* pPtiStart)
447{
448 // OK, we have flagged the start point of this new coastline-normal profile, so create it. Make the start of the profile the centroid of the actual cell that is marked as coast (not the cell under the smoothed vector coast, they may well be different)
449 CGeom2DPoint PtStart; // In external CRS
450 PtStart.SetX(dGridCentroidXToExtCRSX(pPtiStart->nGetX()));
451 PtStart.SetY(dGridCentroidYToExtCRSY(pPtiStart->nGetY()));
452
453 CGeom2DPoint PtEnd; // In external CRS
454 CGeom2DIPoint PtiEnd; // In grid CRS
455 int const nRet = nGetCoastNormalEndPoint(nCoast, nProfileStartPoint, nCoastSize, &PtStart, m_dCoastNormalLength, &PtEnd, &PtiEnd, bIntervention);
456
458 {
459 // Could not solve end-point equation, so forget about this profile
460 return nRet;
461 }
462
463 int const nXEnd = PtiEnd.nGetX();
464 int const nYEnd = PtiEnd.nGetY();
465
466 // Safety check: is the end point in the contiguous sea?
467 if (! m_pRasterGrid->m_Cell[nXEnd][nYEnd].bIsInContiguousSea())
468 {
469 // if (m_nLogFileDetail >= LOG_FILE_ALL)
470 // LogStream << m_ulIter << ": coast " << nCoast << ", possible profile with start point " << nProfileStartPoint << " has inland end point at [" << nXEnd << "][" << nYEnd << "] = {" << dGridCentroidXToExtCRSX(nXEnd) << ", " << dGridCentroidYToExtCRSY(nYEnd) << "}, ignoring" << endl;
471
473 }
474
475 // Safety check: is the water depth at the end point less than the depth of closure?
476 if (m_pRasterGrid->m_Cell[nXEnd][nYEnd].dGetSeaDepth() < m_dDepthOfClosure)
477 {
478 // if (m_nLogFileDetail >= LOG_FILE_ALL)
479 // LogStream << m_ulIter << ": coast " << nCoast << ", possible profile with start point " << nProfileStartPoint << " is too short for depth of closure " << m_dDepthOfClosure << " at end point [" << nXEnd << "][" << nYEnd << "] = {" << dGridCentroidXToExtCRSX(nXEnd) << ", " << dGridCentroidYToExtCRSY(nYEnd) << "}, ignoring" << endl;
480
482 }
483
484 // No problems, so create the new profile
485 CGeomProfile* pProfile = new CGeomProfile(nCoast, nProfileStartPoint, nProfile, pPtiStart, &PtiEnd, bIntervention);
486
487 // And create the profile's coastline-normal vector. Only two points (start and end points, both external CRS) are stored
488 vector<CGeom2DPoint> VNormal;
489 VNormal.push_back(PtStart);
490 VNormal.push_back(PtEnd);
491
492 // Set the start and end points (external CRS) of the profile
493 pProfile->SetPointsInProfile(&VNormal);
494
495 // Create the profile's CGeomMultiLine then set nProfile as the only co-incident profile of the only line segment
496 pProfile->AppendLineSegment();
497 pProfile->AppendCoincidentProfileToLineSegments(make_pair(nProfile, 0));
498
499 // Save the profile, note that several fields in the profile are still blank
500 m_VCoast[nCoast].AppendProfile(pProfile);
501
502 // // DEBUG CODE =================
503 // LogStream << "in nCreateProfile() ===========" << endl;
504 // // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(nProfile);
505 // LogStream << pProfile->nGetProfileID() << "\t";
506 //
507 // int nPointsInProfile = pProfile->nGetProfileSize();
508 //
509 // for (int nPoint = 0; nPoint < nPointsInProfile; nPoint++)
510 // {
511 // CGeom2DPoint Pt = *pProfile->pPtGetPointInProfile(nPoint);
512 // LogStream << " {" << Pt.dGetX() << ", " << Pt.dGetY() << "}";
513 // }
514 // LogStream << endl << "===========" << endl;
515 // // DEBUG CODE =================
516
517 // assert(pProfile->nGetProfileSize() > 0);
518
519 LogStream << m_ulIter << ": coast " << nCoast << " profile " << nProfile << " (nCoastID = " << pProfile->nGetProfileID() << ") created at coast point " << nProfileStartPoint << " from [" << pPtiStart->nGetX() << "][" << pPtiStart->nGetY() << "] = {" << PtStart.dGetX() << ", " << PtStart.dGetY() << "} to [" << PtiEnd.nGetX() << "][" << PtiEnd.nGetY() << "] = {" << PtEnd.dGetX() << ", " << PtEnd.dGetY() << "}" << (pProfile->bIsIntervention() ? ", from intervention" : "") << endl;
520
521 return RTN_OK;
522}
523
524//===============================================================================================================================
526//===============================================================================================================================
527int CSimulation::nLocateAndCreateGridEdgeProfile(bool const bCoastStart, int const nCoast, int& nProfile)
528{
529 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
530 int const nHandedness = m_VCoast[nCoast].nGetSeaHandedness();
531 int const nProfileLen = nRound(m_dCoastNormalLength / m_dCellSide); // Profile length in grid CRS
532 int nProfileStartEdge;
533
534 CGeom2DIPoint PtiProfileStart; // In grid CRS
535 vector<CGeom2DIPoint> VPtiNormalPoints; // In grid CRS
536
537 if (bCoastStart)
538 {
539 // At start of coast
540 PtiProfileStart = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(0); // Grid CRS
541 nProfileStartEdge = m_VCoast[nCoast].nGetStartEdge();
542 }
543 else
544 {
545 // At end of coast
546 PtiProfileStart = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nCoastSize - 1); // Grid CRS
547 nProfileStartEdge = m_VCoast[nCoast].nGetEndEdge();
548 }
549
550 VPtiNormalPoints.push_back(PtiProfileStart);
551
552 // Find the start cell in the list of edge cells
553 auto it = find(m_VEdgeCell.begin(), m_VEdgeCell.end(), PtiProfileStart);
554
555 if (it == m_VEdgeCell.end())
556 {
557 // Not found. This can happen because of rounding problems, i.e. the cell which was stored as the first cell of the raster coastline
559 LogStream << m_ulIter << ": " << ERR << " when constructing start-of-coast profile, [" << PtiProfileStart.nGetX() << "][" << PtiProfileStart.nGetY() << "] = {" << dGridCentroidXToExtCRSX(PtiProfileStart.nGetX()) << ", " << dGridCentroidYToExtCRSY(PtiProfileStart.nGetY()) << "} not found in list of edge cells" << endl;
560
562 }
563
564 // Found
565 int nPos = static_cast<int>(it - m_VEdgeCell.begin());
566
567 // Now construct the edge profile, searching for edge cells
568 for (int n = 0; n < nProfileLen; n++)
569 {
570 if (bCoastStart)
571 {
572 // At start of coast
573 if (nHandedness == LEFT_HANDED)
574 {
575 // The list of edge cells is in clockwise sequence, go in this direction
576 nPos++;
577
578 if (nPos >= static_cast<int>(m_VEdgeCell.size()))
579 {
580 // We've reached the end of the list of edge cells before the profile is long enough. OK, we can live with this
581 break;
582 }
583 }
584 else // Right-handed
585 {
586 // The list of edge cells is in clockwise sequence, go in the opposite direction
587 nPos--;
588
589 if (nPos < 0)
590 {
591 // We've reached the beginning of the list of edge cells before the profile is long enough. OK, we can live with this
592 break;
593 }
594 }
595 }
596 else
597 {
598 // At end of coast
599 if (nHandedness == LEFT_HANDED)
600 {
601 // The list of edge cells is in clockwise sequence, go in the opposite direction
602 nPos--;
603
604 if (nPos < 0)
605 {
606 // We've reached the beginning of the list of edge cells before the profile is long enough. OK, we can live with this
607 break;
608 }
609 }
610 else // Right-handed
611 {
612 // The list of edge cells is in clockwise sequence, go in this direction
613 nPos++;
614
615 if (nPos >= static_cast<int>(m_VEdgeCell.size()))
616 {
617 // We've reached the end of the list of edge cells before the profile is long enough. OK, we can live with this
618 break;
619 }
620 }
621 }
622
623 if (m_VEdgeCellEdge[nPos] != nProfileStartEdge)
624 {
625 // We've reached the end of a grid side before the profile is long enough. OK, we can live with this
626 break;
627 }
628
629 // All OK, so append this grid-edge cell, making sure that there is no gap between this and the previously-appended cell (if there is, will get problems with cell-by-cell fill)
630 AppendEnsureNoGap(&VPtiNormalPoints, &m_VEdgeCell[nPos]);
631 }
632
633 int nProfileStartPoint;
634 CGeomProfile* pProfile;
635 CGeom2DIPoint const PtiDummy(INT_NODATA, INT_NODATA);
636
637 if (bCoastStart)
638 {
639 nProfileStartPoint = 0;
640
641 // Create the new start-of-coast profile
642 pProfile = new CGeomProfile(nCoast, nProfileStartPoint, nProfile, &PtiProfileStart, &PtiDummy, false);
643
644 // Mark this as a start-of-coast profile
645 pProfile->SetStartOfCoast(true);
646 }
647 else
648 {
649 nProfileStartPoint = nCoastSize - 1;
650
651 // Create the new end-of-coast profile
652 pProfile = new CGeomProfile(nCoast, nProfileStartPoint, nProfile, &PtiProfileStart, &PtiDummy, false);
653
654 // Mark this as an end-of-coast profile
655 pProfile->SetEndOfCoast(true);
656 }
657
658 // Create the list of cells 'under' this grid-edge profile. Note that more than two cells are stored
659 for (unsigned int n = 0; n < VPtiNormalPoints.size(); n++)
660 {
661 int const nX = VPtiNormalPoints[n].nGetX();
662 int const nY = VPtiNormalPoints[n].nGetY();
663
664 // Mark each cell in the raster grid
665 m_pRasterGrid->m_Cell[nX][nY].SetCoastAndProfileID(nCoast, nProfile);
666
667 // Store the raster grid coordinates in the profile object
668 pProfile->AppendCellInProfile(nX, nY);
669
670 CGeom2DPoint const Pt(dGridCentroidXToExtCRSX(nX), dGridCentroidYToExtCRSY(nY)); // In external CRS
671
672 // Store the external coordinates in the profile object. Note that for this grid-edge profile, the coordinates of the cells and the coordinates of points on the profile itself are identical, this is not the case for ordinary profiles
673 pProfile->AppendPointInProfile(&Pt);
674 }
675
676 int const nEndX = VPtiNormalPoints.back().nGetX();
677 int const nEndY = VPtiNormalPoints.back().nGetY();
678 CGeom2DIPoint const PtiEnd(nEndX, nEndY);
679 pProfile->SetEndPoint(&PtiEnd);
680
681 // Get the deep water wave height and orientation values at the end of the profile
682 double const dDeepWaterWaveHeight = m_pRasterGrid->m_Cell[nEndX][nEndY].dGetCellDeepWaterWaveHeight();
683 double const dDeepWaterWaveAngle = m_pRasterGrid->m_Cell[nEndX][nEndY].dGetCellDeepWaterWaveAngle();
684 double const dDeepWaterWavePeriod = m_pRasterGrid->m_Cell[nEndX][nEndY].dGetCellDeepWaterWavePeriod();
685
686 // And store them in this profile
687 pProfile->SetProfileDeepWaterWaveHeight(dDeepWaterWaveHeight);
688 pProfile->SetProfileDeepWaterWaveAngle(dDeepWaterWaveAngle);
689 pProfile->SetProfileDeepWaterWavePeriod(dDeepWaterWavePeriod);
690
691 // Create the profile's CGeomMultiLine then set nProfile as the only co-incident profile of the only line segment
692 pProfile->AppendLineSegment();
693 pProfile->AppendCoincidentProfileToLineSegments(make_pair(nProfile, 0));
694
695 // Store the grid-edge profile
696 m_VCoast[nCoast].AppendProfile(pProfile);
697 m_VCoast[nCoast].SetProfileAtCoastPoint(nProfileStartPoint, pProfile);
698
700 LogStream << m_ulIter << ": coast " << nCoast << " grid-edge profile " << nProfile << " (nCoastID = " << pProfile->nGetProfileID() << ") created at coast " << (bCoastStart ? "start" : "end") << " point " << (bCoastStart ? 0 : nCoastSize - 1) << ", from [" << PtiProfileStart.nGetX() << "][" << PtiProfileStart.nGetY() << "] = {" << dGridCentroidXToExtCRSX(PtiProfileStart.nGetX()) << ", " << dGridCentroidYToExtCRSY(PtiProfileStart.nGetY()) << "} to [" << VPtiNormalPoints.back().nGetX() << "][" << VPtiNormalPoints.back().nGetY() << "] = {" << dGridCentroidXToExtCRSX(VPtiNormalPoints.back().nGetX()) << ", " << dGridCentroidYToExtCRSY(VPtiNormalPoints.back().nGetY()) << "} profile length = " << VPtiNormalPoints.size() << endl;
701
702 // assert(pProfile->nGetProfileSize() > 0);
703
704 return RTN_OK;
705}
706
707//===============================================================================================================================
709//===============================================================================================================================
710int CSimulation::nGetCoastNormalEndPoint(int const nCoast, int const nStartCoastPoint, int const nCoastSize, CGeom2DPoint const* pPtStart, double const dLineLength, CGeom2DPoint* pPtEnd, CGeom2DIPoint* pPtiEnd, bool const bIntervention)
711{
712 int const AVGSIZE = 21; // TODO 011 This should be a user input
713
714 CGeom2DPoint PtBefore;
715 CGeom2DPoint PtAfter;
716
717 if (bIntervention)
718 {
719 // This is an intervention profile, so just use one point on either side (coordinates in external CRS). Note this this assumes that this intervention profile is not at the start or end of the coastline
720 PtBefore = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nStartCoastPoint - 1);
721 PtAfter = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nStartCoastPoint + 1);
722 }
723
724 else
725 {
726 // This is not an intervention, so put a maximum of AVGSIZE points before the start point into a vector
727 vector<CGeom2DPoint> PtBeforeToAverage;
728
729 for (int n = 1; n <= AVGSIZE; n++)
730 {
731 int const nPoint = nStartCoastPoint - n;
732
733 if (nPoint < 0)
734 break;
735
736 PtBeforeToAverage.push_back(*m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nPoint));
737 }
738
739 // Put a maximum of AVGSIZE points after the start point into a vector
740 vector<CGeom2DPoint> PtAfterToAverage;
741
742 for (int n = 1; n <= AVGSIZE; n++)
743 {
744 int const nPoint = nStartCoastPoint + n;
745
746 if (nPoint > nCoastSize - 1)
747 break;
748
749 PtAfterToAverage.push_back(*m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nPoint));
750 }
751
752 // Now average each of these vectors of points: results are in PtBefore and PtAfter (coordinates in external CRS)
753 PtBefore = PtAverage(&PtBeforeToAverage);
754 PtAfter = PtAverage(&PtAfterToAverage);
755 }
756
757 // Get the y = a * x + b equation of the straight line linking the coastline points before and after 'this' coastline point. For this linking line, slope a = (y2 - y1) / (x2 - x1)
758 double const dYDiff = PtAfter.dGetY() - PtBefore.dGetY();
759 double const dXDiff = PtAfter.dGetX() - PtBefore.dGetX();
760
761 double dXEnd1 = 0, dXEnd2 = 0, dYEnd1 = 0, dYEnd2 = 0;
762
763 if (bFPIsEqual(dYDiff, 0.0, TOLERANCE))
764 {
765 // The linking line runs W-E or E-W, so a straight line at right angles to this runs N-S or S-N. Calculate the two possible end points for this coastline-normal profile
766 dXEnd1 = dXEnd2 = pPtStart->dGetX();
767 dYEnd1 = pPtStart->dGetY() + dLineLength;
768 dYEnd2 = pPtStart->dGetY() - dLineLength;
769 }
770
771 else if (bFPIsEqual(dXDiff, 0.0, TOLERANCE))
772 {
773 // The linking line runs N-S or S-N, so a straight line at right angles to this runs W-E or E-W. Calculate the two possible end points for this coastline-normal profile
774 dYEnd1 = dYEnd2 = pPtStart->dGetY();
775 dXEnd1 = pPtStart->dGetX() + dLineLength;
776 dXEnd2 = pPtStart->dGetX() - dLineLength;
777 }
778
779 else
780 {
781 // The linking line runs neither W-E nor N-S so we have to work a bit harder to find the end-point of the coastline-normal profile
782 double const dA = dYDiff / dXDiff;
783
784 // Now calculate the equation of the straight line which is perpendicular to this linking line
785 double const dAPerp = -1 / dA;
786 double const dBPerp = pPtStart->dGetY() - (dAPerp * pPtStart->dGetX());
787
788 // Calculate the end point of the profile: first do some substitution then rearrange as a quadratic equation i.e. in the form Ax^2 + Bx + C = 0 (see http://math.stackexchange.com/questions/228841/how-do-i-calculate-the-intersections-of-a-straight-line-and-a-circle)
789 double const dQuadA = 1 + (dAPerp * dAPerp);
790 double const dQuadB = 2 * ((dBPerp * dAPerp) - (dAPerp * pPtStart->dGetY()) - pPtStart->dGetX());
791 double const dQuadC = ((pPtStart->dGetX() * pPtStart->dGetX()) + (pPtStart->dGetY() * pPtStart->dGetY()) + (dBPerp * dBPerp) - (2 * pPtStart->dGetY() * dBPerp) - (dLineLength * dLineLength));
792
793 // Solve for x and y using the quadratic formula x = (−B ± sqrt(B^2 − 4AC)) / 2A
794 double const dDiscriminant = (dQuadB * dQuadB) - (4 * dQuadA * dQuadC);
795
796 if (dDiscriminant < 0)
797 {
798 LogStream << ERR << "timestep " << m_ulIter << ": discriminant < 0 when finding profile end point on coastline " << nCoast << ", from coastline point " << nStartCoastPoint << "), ignored" << endl;
800 }
801
802 dXEnd1 = (-dQuadB + sqrt(dDiscriminant)) / (2 * dQuadA);
803 dYEnd1 = (dAPerp * dXEnd1) + dBPerp;
804 dXEnd2 = (-dQuadB - sqrt(dDiscriminant)) / (2 * dQuadA);
805 dYEnd2 = (dAPerp * dXEnd2) + dBPerp;
806 }
807
808 // We have two possible solutions, so decide which of the two endpoints to use then create the profile end-point (coordinates in external CRS)
809 int const nSeaHand = m_VCoast[nCoast].nGetSeaHandedness(); // Assumes handedness is either 0 or 1 (i.e. not -1)
810 *pPtEnd = PtChooseEndPoint(nSeaHand, &PtBefore, &PtAfter, dXEnd1, dYEnd1, dXEnd2, dYEnd2);
811
812 // Check that pPtiEnd is not off the grid. Note that pPtiEnd is not necessarily a cell centroid
813 pPtiEnd->SetXY(nRound(dExtCRSXToGridX(pPtEnd->dGetX())), nRound(dExtCRSYToGridY(pPtEnd->dGetY())));
814
815 if (! bIsWithinValidGrid(pPtiEnd))
816 {
817 // LogStream << m_ulIter << ": profile endpoint is outside grid [" << pPtiEnd->nGetX() << "][" << pPtiEnd->nGetY() << "] = {" << pPtEnd->dGetX() << ", " << pPtEnd->dGetY() << "}. The profile starts at coastline point " << nStartCoastPoint << " = {" << pPtStart->dGetX() << ", " << pPtStart->dGetY() << "}" << endl;
818
819 // The end point is off the grid, so constrain it to be within the valid grid
820 CGeom2DIPoint const PtiStart(nRound(dExtCRSXToGridX(pPtStart->dGetX())), nRound(dExtCRSYToGridY(pPtStart->dGetY())));
821 KeepWithinValidGrid(&PtiStart, pPtiEnd);
822
823 pPtEnd->SetX(dGridCentroidXToExtCRSX(pPtiEnd->nGetX()));
824 pPtEnd->SetY(dGridCentroidYToExtCRSY(pPtiEnd->nGetY()));
825
826 // LogStream << m_ulIter << ": profile endpoint is now within the grid [" << pPtiEnd->nGetX() << "][" << pPtiEnd->nGetY() << "] = {" << pPtEnd->dGetX() << ", " << pPtEnd->dGetY() << "}. The profile starts at coastline point " << nStartCoastPoint << " = {" << pPtStart->dGetX() << ", " << pPtStart->dGetY() << "}" << endl;
827
828 // return RTN_ERR_PROFILE_ENDPOINT_AT_GRID_EDGE;
829 }
830
831 return RTN_OK;
832}
833
834//===============================================================================================================================
836//===============================================================================================================================
837CGeom2DPoint CSimulation::PtChooseEndPoint(int const nHand, CGeom2DPoint const* PtBefore, CGeom2DPoint const* PtAfter, double const dXEnd1, double const dYEnd1, double const dXEnd2, double const dYEnd2)
838{
839 CGeom2DPoint PtChosen;
840
841 // All coordinates here are in the external CRS, so the origin of the grid is the bottom left
842 if (nHand == RIGHT_HANDED)
843 {
844 // The sea is to the right of the linking line. So which way is the linking line oriented? First check the N-S component
845 if (PtAfter->dGetY() > PtBefore->dGetY())
846 {
847 // We are going S to N and the sea is to the right: the normal endpoint is to the E. We want the larger of the two x values
848 if (dXEnd1 > dXEnd2)
849 {
850 PtChosen.SetX(dXEnd1);
851 PtChosen.SetY(dYEnd1);
852 }
853
854 else
855 {
856 PtChosen.SetX(dXEnd2);
857 PtChosen.SetY(dYEnd2);
858 }
859 }
860
861 else if (PtAfter->dGetY() < PtBefore->dGetY())
862 {
863 // We are going N to S and the sea is to the right: the normal endpoint is to the W. We want the smaller of the two x values
864 if (dXEnd1 < dXEnd2)
865 {
866 PtChosen.SetX(dXEnd1);
867 PtChosen.SetY(dYEnd1);
868 }
869
870 else
871 {
872 PtChosen.SetX(dXEnd2);
873 PtChosen.SetY(dYEnd2);
874 }
875 }
876
877 else
878 {
879 // No N-S component i.e. the linking line is exactly W-E. So check the W-E component
880 if (PtAfter->dGetX() > PtBefore->dGetX())
881 {
882 // We are going W to E and the sea is to the right: the normal endpoint is to the s. We want the smaller of the two y values
883 if (dYEnd1 < dYEnd2)
884 {
885 PtChosen.SetX(dXEnd1);
886 PtChosen.SetY(dYEnd1);
887 }
888
889 else
890 {
891 PtChosen.SetX(dXEnd2);
892 PtChosen.SetY(dYEnd2);
893 }
894 }
895
896 else // Do not check for (PtAfter->dGetX() == PtBefore->dGetX()), since this would mean the two points are co-incident
897 {
898 // We are going E to W and the sea is to the right: the normal endpoint is to the N. We want the larger of the two y values
899 if (dYEnd1 > dYEnd2)
900 {
901 PtChosen.SetX(dXEnd1);
902 PtChosen.SetY(dYEnd1);
903 }
904
905 else
906 {
907 PtChosen.SetX(dXEnd2);
908 PtChosen.SetY(dYEnd2);
909 }
910 }
911 }
912 }
913
914 else // nHand == LEFT_HANDED
915 {
916 // The sea is to the left of the linking line. So which way is the linking line oriented? First check the N-S component
917 if (PtAfter->dGetY() > PtBefore->dGetY())
918 {
919 // We are going S to N and the sea is to the left: the normal endpoint is to the W. We want the smaller of the two x values
920 if (dXEnd1 < dXEnd2)
921 {
922 PtChosen.SetX(dXEnd1);
923 PtChosen.SetY(dYEnd1);
924 }
925
926 else
927 {
928 PtChosen.SetX(dXEnd2);
929 PtChosen.SetY(dYEnd2);
930 }
931 }
932
933 else if (PtAfter->dGetY() < PtBefore->dGetY())
934 {
935 // We are going N to S and the sea is to the left: the normal endpoint is to the E. We want the larger of the two x values
936 if (dXEnd1 > dXEnd2)
937 {
938 PtChosen.SetX(dXEnd1);
939 PtChosen.SetY(dYEnd1);
940 }
941
942 else
943 {
944 PtChosen.SetX(dXEnd2);
945 PtChosen.SetY(dYEnd2);
946 }
947 }
948
949 else
950 {
951 // No N-S component i.e. the linking line is exactly W-E. So check the W-E component
952 if (PtAfter->dGetX() > PtBefore->dGetX())
953 {
954 // We are going W to E and the sea is to the left: the normal endpoint is to the N. We want the larger of the two y values
955 if (dYEnd1 > dYEnd2)
956 {
957 PtChosen.SetX(dXEnd1);
958 PtChosen.SetY(dYEnd1);
959 }
960
961 else
962 {
963 PtChosen.SetX(dXEnd2);
964 PtChosen.SetY(dYEnd2);
965 }
966 }
967
968 else // Do not check for (PtAfter->dGetX() == PtBefore->dGetX()), since this would mean the two points are co-incident
969 {
970 // We are going E to W and the sea is to the left: the normal endpoint is to the S. We want the smaller of the two y values
971 if (dYEnd1 < dYEnd2)
972 {
973 PtChosen.SetX(dXEnd1);
974 PtChosen.SetY(dYEnd1);
975 }
976
977 else
978 {
979 PtChosen.SetX(dXEnd2);
980 PtChosen.SetY(dYEnd2);
981 }
982 }
983 }
984 }
985
986 return PtChosen;
987}
988
989//===============================================================================================================================
991//===============================================================================================================================
993{
994 // Do once for every coastline object
995 int const nCoastLines = static_cast<int>(m_VCoast.size());
996
997 for (int nCoast = 0; nCoast < nCoastLines; nCoast++)
998 {
999 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
1000
1001 // Do once for every profile, in along-coast sequence
1002 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
1003 {
1004 if (! m_VCoast[nCoast].bIsProfileAtCoastPoint(nCoastPoint))
1005 continue;
1006
1007 // There is a profile at this coast point
1008 CGeomProfile* pFirstProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(nCoastPoint);
1009 int const nFirstProfile = pFirstProfile->nGetProfileID();
1010
1011 // Only check this profile if it is problem free, and is not a start- or end-of-coast profile. Continue checking if it has been truncated, however
1012 if (! pFirstProfile->bProfileOKIncTruncated())
1013 {
1014 // LogStream << m_ulIter << ": nCoastPoint = " << nCoastPoint << " pFirstProfile = " << pFirstProfile->nGetProfileID() << " is not OK (could be a start- or end-of-coast profile), abandoning" << endl;
1015 continue;
1016 }
1017
1018 // OK we have found a first profile. Now go along the coast in alternate directions: first down-coast (in the direction of increasing coast point numbers) then up-coast
1019 for (int nDirection = DIRECTION_DOWNCOAST; nDirection <= DIRECTION_UPCOAST; nDirection++)
1020 {
1021 int nStartPoint;
1022
1023 if (nDirection == DIRECTION_DOWNCOAST)
1024 nStartPoint = nCoastPoint + 1;
1025 else
1026 nStartPoint = nCoastPoint - 1;
1027
1028 for (int nSecondCoastPoint = nStartPoint; (nDirection == DIRECTION_DOWNCOAST) ? nSecondCoastPoint < nCoastSize : nSecondCoastPoint >= 0; (nDirection == DIRECTION_DOWNCOAST) ? nSecondCoastPoint++ : nSecondCoastPoint--)
1029 // // In this direction, look at profiles which are increasingly close to the first profile
1030 // int nStartPoint;
1031 // if (nDirection == DIRECTION_DOWNCOAST)
1032 // nStartPoint = 0;
1033 // else
1034 // nStartPoint = nCoastSize - 1;
1035 //
1036 // for (int nSecondCoastPoint = nStartPoint; (nDirection == DIRECTION_DOWNCOAST) ? nSecondCoastPoint < nCoastPoint : nSecondCoastPoint > nCoastPoint; (nDirection == DIRECTION_DOWNCOAST) ? nSecondCoastPoint++ : nSecondCoastPoint--)
1037 {
1038 if (m_VCoast[nCoast].bIsProfileAtCoastPoint(nSecondCoastPoint))
1039 {
1040 // There is a profile at the second coast point, so get a pointer to it
1041 CGeomProfile* pSecondProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(nSecondCoastPoint);
1042 int const nSecondProfile = pSecondProfile->nGetProfileID();
1043
1044 // LogStream << m_ulIter << ": " << (nDirection == DIRECTION_DOWNCOAST ? "down" : "up") << "-coast search, nCoastPoint = " << nCoastPoint << " nSecondCoastPoint = " << nSecondCoastPoint << " (profiles " << pFirstProfile->nGetProfileID() << " and " << pSecondProfile->nGetProfileID() << ")" << endl;
1045
1046 // Only check this profile if it is problem free, and is not a start- or end-of-coast profile. Continue checking if it has been truncated, however
1047 if (! pSecondProfile->bProfileOKIncTruncated())
1048 {
1049 // LogStream << m_ulIter << ": second profile = " << pSecondProfile->nGetProfileID() << " is not OK (could be a start- or end-of-coast profile), abandoning" << endl;
1050 continue;
1051 }
1052
1053 // Only check these two profiles for intersection if they are are not co-incident in the final line segment of both profiles (i.e. the profiles have not already intersected)
1054 if ((pFirstProfile->bFindProfileInCoincidentProfilesOfLastLineSegment(nSecondProfile)) || (pSecondProfile->bFindProfileInCoincidentProfilesOfLastLineSegment(nFirstProfile)))
1055 {
1056 // LogStream << m_ulIter << ": profiles " << pFirstProfile->nGetProfileID() << " and " << pSecondProfile->nGetProfileID() << " are are not co-incident in the final line segment of both profiles (i.e. the profiles have not already intersected), abandoning" << endl;
1057 continue;
1058 }
1059
1060 // OK go for it
1061 int nProf1LineSeg = 0;
1062 int nProf2LineSeg = 0;
1063 double dIntersectX = 0;
1064 double dIntersectY = 0;
1065 double dAvgEndX = 0;
1066 double dAvgEndY = 0;
1067
1068 if (bCheckForIntersection(pFirstProfile, pSecondProfile, nProf1LineSeg, nProf2LineSeg, dIntersectX, dIntersectY, dAvgEndX, dAvgEndY))
1069 {
1070 // The profiles intersect. Decide which profile to truncate, and which to retain
1071 int nPoint = -1;
1072
1073 if (pFirstProfile->bIsIntervention())
1074 {
1075 LogStream << m_ulIter << ": profiles " << nFirstProfile << " and " << nSecondProfile << " intersect, truncate " << nFirstProfile << " since it is an intervention profile" << endl;
1076
1077 // Truncate the first profile, since it is an intervention profile
1078 TruncateOneProfileRetainOtherProfile(nCoast, pFirstProfile, pSecondProfile, dIntersectX, dIntersectY, nProf1LineSeg, nProf2LineSeg, false);
1079 }
1080
1081 else if (pSecondProfile->bIsIntervention())
1082 {
1083 LogStream << m_ulIter << ": profiles " << nFirstProfile << " and " << nSecondProfile << " intersect, truncate " << nSecondProfile << " since it is an intervention profile" << endl;
1084
1085 // Truncate the second profile, since it is an intervention profile
1086 TruncateOneProfileRetainOtherProfile(nCoast, pSecondProfile, pFirstProfile, dIntersectX, dIntersectY, nProf2LineSeg, nProf1LineSeg, false);
1087 }
1088
1089 // Is the point of intersection already present in the first profile (i.e. because there has already been an intersection at this point between the first profile and some other profile)?
1090 else if (pFirstProfile->bIsPointInProfile(dIntersectX, dIntersectY, nPoint))
1091 {
1092 LogStream << m_ulIter << ": profiles " << nFirstProfile << " and " << nSecondProfile << " intersect, but point {" << dIntersectX << ", " << dIntersectY << "} is already present in profile " << nFirstProfile << " as point " << nPoint << endl;
1093
1094 // Truncate the second profile and merge it with the first profile
1095 TruncateOneProfileRetainOtherProfile(nCoast, pSecondProfile, pFirstProfile, dIntersectX, dIntersectY, nProf2LineSeg, nProf1LineSeg, true);
1096 }
1097
1098 // Is the point of intersection already present in the second profile?
1099 else if (pSecondProfile->bIsPointInProfile(dIntersectX, dIntersectY, nPoint))
1100 {
1101 LogStream << m_ulIter << ": profiles " << nFirstProfile << " and " << nSecondProfile << " intersect, but point {" << dIntersectX << ", " << dIntersectY << "} is already present in profile " << nSecondProfile << " as point " << nPoint << endl;
1102
1103 // Truncate the first profile and merge it with the second profile
1104 TruncateOneProfileRetainOtherProfile(nCoast, pFirstProfile, pSecondProfile, dIntersectX, dIntersectY, nProf1LineSeg, nProf2LineSeg, true);
1105 }
1106
1107 else
1108 {
1109 // The point of intersection is not already present in either profile, so get the number of line segments of each profile
1110 int const nFirstProfileLineSegments = pFirstProfile->nGetNumLineSegments();
1111 int const nSecondProfileLineSegments = pSecondProfile->nGetNumLineSegments();
1112
1113 // assert(nProf1LineSeg < nFirstProfileLineSegments);
1114 // assert(nProf2LineSeg < nSecondProfileLineSegments);
1115
1116 // Next check whether the point of intersection is on the final line segment of both profiles
1117 if ((nProf1LineSeg == (nFirstProfileLineSegments - 1)) && (nProf2LineSeg == (nSecondProfileLineSegments - 1)))
1118 {
1119 // Yes, the point of intersection is on the final line segment of both profiles, so merge the profiles seaward of the point of intersection
1120 MergeProfilesAtFinalLineSegments(nCoast, pFirstProfile, pSecondProfile, nFirstProfileLineSegments, nSecondProfileLineSegments, dIntersectX, dIntersectY, dAvgEndX, dAvgEndY);
1121
1122 // LogStream << m_ulIter << ": " << ((nDirection == DIRECTION_DOWNCOAST) ? "down" : "up") << "-coast search, end-segment intersection between profiles " << nFirstProfile << " and " << nSecondProfile << " at [" << dIntersectX << ", " << dIntersectY << "] in line segment [" << nProf1LineSeg << "] of " << nFirstProfileLineSegments << " segments, and line segment [" << nProf2LineSeg << "] of " << nSecondProfileLineSegments << " segments, respectively" << endl;
1123
1124 // // DEBUG CODE =============================================================================================
1125 // int nSizeTmp = pFirstProfile->nGetProfileSize();
1126 // CGeom2DPoint PtEndTmp = *pFirstProfile->pPtGetPointInProfile(nSizeTmp-1);
1127 //
1128 // LogStream << m_ulIter << ": end of first profile (" << nFirstProfile << ") is point " << nSizeTmp-1 << " at [" << dExtCRSXToGridX(PtEndTmp.dGetX()) << "][" << dExtCRSYToGridY(PtEndTmp.dGetY()) << "} = {" << PtEndTmp.dGetX() << ", " << PtEndTmp.dGetY() << "}" << endl;
1129 //
1130 // nSizeTmp = pSecondProfile->nGetProfileSize();
1131 // PtEndTmp = *pSecondProfile->pPtGetPointInProfile(nSizeTmp-1);
1132 //
1133 // LogStream << m_ulIter << ": end of second profile (" << nSecondProfile << ") is point " << nSizeTmp-1 << " at [" << dExtCRSXToGridX(PtEndTmp.dGetX()) << "][" << dExtCRSYToGridY(PtEndTmp.dGetY()) << "} = {" << PtEndTmp.dGetX() << ", " << PtEndTmp.dGetY() << "}" << endl;
1134 // // DEBUG CODE =============================================================================================
1135 }
1136
1137 else
1138 {
1139 // The profiles intersect, but the point of intersection is not on the final line segment of both profiles. One of the profiles will be truncated, the other profile will be retained
1140 // LogStream << m_ulIter << ": " << ((nDirection == DIRECTION_DOWNCOAST) ? "down" : "up") << "-coast search, intersection (NOT both end segments) between profiles " << nFirstProfile << " and " << nSecondProfile << " at [" << dIntersectX << ", " << dIntersectY << "] in line segment [" << nProf1LineSeg << "] of " << nFirstProfileLineSegments << ", and line segment [" << nProf2LineSeg << "] of " << nSecondProfileLineSegments << ", respectively" << endl;
1141
1142 // Decide which profile to truncate, and which to retain
1143 if (pFirstProfile->bIsIntervention())
1144 {
1145 // Truncate the first profile, since it is an intervention profile
1146 // LogStream << m_ulIter << ": pFirstProfile = " << pFirstProfile->nGetProfileID() << " pSecondProfile = " << pSecondProfile->nGetProfileID() << ", pFirstProfile is an intervention profile, so truncate pFirstProfile (" << pFirstProfile->nGetProfileID() << ")" << endl;
1147
1148 TruncateOneProfileRetainOtherProfile(nCoast, pFirstProfile, pSecondProfile, dIntersectX, dIntersectY, nProf1LineSeg, nProf2LineSeg, false);
1149 }
1150
1151 else if (pSecondProfile->bIsIntervention())
1152 {
1153 // Truncate the second profile, since it is an intervention profile
1154 // LogStream << m_ulIter << ": pFirstProfile = " << pFirstProfile->nGetProfileID() << " pSecondProfile = " << pSecondProfile->nGetProfileID() << ", pSecondProfile is an intervention profile, so truncate pSecondProfile (" << pSecondProfile->nGetProfileID() << ")" << endl;
1155
1156 TruncateOneProfileRetainOtherProfile(nCoast, pSecondProfile, pFirstProfile, dIntersectX, dIntersectY, nProf2LineSeg, nProf1LineSeg, false);
1157 }
1158
1159 else if (nFirstProfileLineSegments < nSecondProfileLineSegments)
1160 {
1161 // Truncate the first profile, since it has a smaller number of line segments
1162 // LogStream << m_ulIter << ": pFirstProfile = " << pFirstProfile->nGetProfileID() << " pSecondProfile = " << pSecondProfile->nGetProfileID() << ", pFirstProfile has a smaller number of line segments, so truncate pFirstProfile (" << pFirstProfile->nGetProfileID() << ")" << endl;
1163
1164 TruncateOneProfileRetainOtherProfile(nCoast, pFirstProfile, pSecondProfile, dIntersectX, dIntersectY, nProf1LineSeg, nProf2LineSeg, false);
1165 }
1166
1167 else if (nFirstProfileLineSegments > nSecondProfileLineSegments)
1168 {
1169 // Truncate the second profile, since it has a smaller number of line segments
1170 // LogStream << m_ulIter << ": pFirstProfile = " << pFirstProfile->nGetProfileID() << " pSecondProfile = " << pSecondProfile->nGetProfileID() << ", pSecondProfile has a smaller number of line segments, so truncate pSecondProfile (" << pSecondProfile->nGetProfileID() << ")" << endl;
1171
1172 TruncateOneProfileRetainOtherProfile(nCoast, pSecondProfile, pFirstProfile, dIntersectX, dIntersectY, nProf2LineSeg, nProf1LineSeg, false);
1173 }
1174
1175 else
1176 {
1177 // Both profiles have the same number of line segments, so choose randomly. Draw a sample from the unit normal distribution using random number generator 1
1178 double const dRand = m_dGetFromUnitNormalDist(m_Rand[0]);
1179
1180 if (dRand >= 0.0)
1181 {
1182 // LogStream << m_ulIter << ": pFirstProfile = " << pFirstProfile->nGetProfileID() << " pSecondProfile = " << pSecondProfile->nGetProfileID() << ", same number of line segment, randomly truncate pFirstProfile" << endl;
1183
1184 TruncateOneProfileRetainOtherProfile(nCoast, pFirstProfile, pSecondProfile, dIntersectX, dIntersectY, nProf1LineSeg, nProf2LineSeg, false);
1185 }
1186
1187 else
1188 {
1189 // LogStream << m_ulIter << ": pFirstProfile = " << pFirstProfile->nGetProfileID() << " pSecondProfile = " << pSecondProfile->nGetProfileID() << ", same number of line segment, randomly truncate pSecondProfile" << endl;
1190
1191 TruncateOneProfileRetainOtherProfile(nCoast, pSecondProfile, pFirstProfile, dIntersectX, dIntersectY, nProf2LineSeg, nProf1LineSeg, false);
1192 }
1193 }
1194 }
1195 }
1196 }
1197 }
1198 }
1199 }
1200 }
1201 }
1202}
1203
1204//===============================================================================================================================
1206//===============================================================================================================================
1208{
1209 // Check to see which coastline-normal profiles intersect. Then modify intersecting profiles so that the sections of each profile seaward of the point of intersection are 'shared' i.e. are multi-lines. This creates the boundaries of the triangular polygons
1211
1212 // Again check the normal profiles for insufficient length: is the water depth at the end point less than the depth of closure? We do this again because some profiles may have been shortened as a result of intersection. Do once for every coastline object
1213 for (unsigned int nCoast = 0; nCoast < m_VCoast.size(); nCoast++)
1214 {
1215 for (int n = 0; n < m_VCoast[nCoast].nGetNumProfiles(); n++)
1216 {
1217 CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(n);
1218 int const nProfile = pProfile->nGetProfileID();
1219
1220 if (pProfile->bProfileOK())
1221 {
1222 int const nSize = pProfile->nGetProfileSize();
1223
1224 // Safety check
1225 if (nSize == 0)
1226 {
1227 // pProfile->SetTooShort(true);
1228 m_VCoast[nCoast].pGetProfile(nProfile)->SetTooShort(true);
1229 LogStream << "Profile " << nProfile << " is too short, size = " << nSize << endl;
1230 continue;
1231 }
1232
1233 CGeom2DPoint const* pPtEnd = pProfile->pPtGetPointInProfile(nSize - 1);
1234 CGeom2DIPoint const PtiEnd = PtiExtCRSToGridRound(pPtEnd);
1235 int nXEnd = PtiEnd.nGetX();
1236 int nYEnd = PtiEnd.nGetY();
1237
1238 // Safety checks: the point may be outside the grid, so keep it within the grid
1239 nXEnd = tMin(nXEnd, m_nXGridSize - 1);
1240 nYEnd = tMin(nYEnd, m_nYGridSize - 1);
1241 nXEnd = tMax(nXEnd, 0);
1242 nYEnd = tMax(nYEnd, 0);
1243
1244 // pProfile->SetEndPoint(&PtiEnd);
1245 m_VCoast[nCoast].pGetProfile(nProfile)->SetEndPoint(&PtiEnd);
1246
1247 if (m_pRasterGrid->m_Cell[nXEnd][nYEnd].dGetSeaDepth() < m_dDepthOfClosure)
1248 {
1250 LogStream << m_ulIter << ": coast " << nCoast << ", profile " << nProfile << " is invalid, is too short for depth of closure " << m_dDepthOfClosure << " at end point [" << nXEnd << "][" << nYEnd << "] = {" << pPtEnd->dGetX() << ", " << pPtEnd->dGetY() << "}, flagging as too short" << endl;
1251
1252 // pProfile->SetTooShort(true);
1253 m_VCoast[nCoast].pGetProfile(nProfile)->SetTooShort(true);
1254 }
1255 }
1256 }
1257
1258 // For this coast, put all valid coastline-normal profiles (apart from the profiles at the start and end of the coast, since they have already been done) onto the raster grid. But if the profile is not long enough, crosses a coastline, hits dry land, or hits another profile, then mark the profile as invalid
1259 int nValidProfiles = 0;
1260 MarkProfilesOnGrid(nCoast, nValidProfiles);
1261
1262 if (nValidProfiles == 0)
1263 {
1264 // Problem! No valid profiles, so quit
1265 cerr << m_ulIter << ": " << ERR << "no coastline-normal profiles created" << endl;
1266 return RTN_ERR_NO_PROFILES_2;
1267 }
1268
1269 // // DEBUG CODE ===========================================================================================================
1270 // if (m_ulIter == 109)
1271 // {
1272 // string strOutFile = m_strOutPath;
1273 // strOutFile += "00_profile_raster_";
1274 // strOutFile += to_string(m_ulIter);
1275 // strOutFile += ".tif";
1276 //
1277 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
1278 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
1279 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
1280 // pDataSet->SetGeoTransform(m_dGeoTransform);
1281 //
1282 // int nn = 0;
1283 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
1284 // for (int nY = 0; nY < m_nYGridSize; nY++)
1285 // {
1286 // for (int nX = 0; nX < m_nXGridSize; nX++)
1287 // {
1288 // if (m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
1289 // pdRaster[nn] = -1;
1290 // else
1291 // {
1292 //
1293 // // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
1294 // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].nGetProfileID();
1295 // }
1296 //
1297 // nn++;
1298 // }
1299 // }
1300 //
1301 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
1302 // pBand->SetNoDataValue(m_dMissingValue);
1303 // int nRet1 = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
1304 //
1305 // if (nRet1 == CE_Failure)
1306 // return RTN_ERR_GRIDCREATE;
1307 //
1308 // GDALClose(pDataSet);
1309 // delete[] pdRaster;
1310 // }
1311 // // DEBUG CODE ===========================================================================================================
1312 }
1313
1314 return RTN_OK;
1315}
1316
1317//===============================================================================================================================
1319//===============================================================================================================================
1320bool CSimulation::bCheckForIntersection(CGeomProfile* const pVProfile1, CGeomProfile* const pVProfile2, int& nProfile1LineSegment, int& nProfile2LineSegment, double& dXIntersect, double& dYIntersect, double& dXAvgEnd, double& dYAvgEnd)
1321{
1322 // For both profiles, look at all line segments
1323 int const nProfile1NumSegments = pVProfile1->nGetNumLineSegments();
1324 int const nProfile2NumSegments = pVProfile2->nGetNumLineSegments();
1325 // nProfile1Size = pVProfile1->nGetProfileSize(),
1326 // nProfile2Size = pVProfile2->nGetProfileSize();
1327
1328 // assert(nProfile1Size == nProfile1NumSegments+1);
1329 // assert(nProfile2Size == nProfile2NumSegments+1);
1330
1331 for (int i = 0; i < nProfile1NumSegments; i++)
1332 {
1333 for (int j = 0; j < nProfile2NumSegments; j++)
1334 {
1335 // In external coordinates
1336 double const dX1 = pVProfile1->pPtVGetPoints()->at(i).dGetX();
1337 double const dY1 = pVProfile1->pPtVGetPoints()->at(i).dGetY();
1338 double const dX2 = pVProfile1->pPtVGetPoints()->at(i + 1).dGetX();
1339 double const dY2 = pVProfile1->pPtVGetPoints()->at(i + 1).dGetY();
1340
1341 double const dX3 = pVProfile2->pPtVGetPoints()->at(j).dGetX();
1342 double const dY3 = pVProfile2->pPtVGetPoints()->at(j).dGetY();
1343 double const dX4 = pVProfile2->pPtVGetPoints()->at(j + 1).dGetX();
1344 double const dY4 = pVProfile2->pPtVGetPoints()->at(j + 1).dGetY();
1345
1346 // Uses Cramer's Rule to solve the equations. Modified from code at http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect (in turn based on Andre LeMothe's "Tricks of the Windows Game Programming Gurus")
1347 double const dDiffX1 = dX2 - dX1;
1348 double const dDiffY1 = dY2 - dY1;
1349 double const dDiffX2 = dX4 - dX3;
1350 double const dDiffY2 = dY4 - dY3;
1351
1352 double dS = -999;
1353 double dT = -999;
1354 double dTmp = 0;
1355
1356 dTmp = -dDiffX2 * dDiffY1 + dDiffX1 * dDiffY2;
1357
1358 if (! bFPIsEqual(dTmp, 0.0, TOLERANCE))
1359 dS = (-dDiffY1 * (dX1 - dX3) + dDiffX1 * (dY1 - dY3)) / dTmp;
1360
1361 dTmp = -dDiffX2 * dDiffY1 + dDiffX1 * dDiffY2;
1362
1363 if (! bFPIsEqual(dTmp, 0.0, TOLERANCE))
1364 dT = (dDiffX2 * (dY1 - dY3) - dDiffY2 * (dX1 - dX3)) / dTmp;
1365
1366 if (dS >= 0 && dS <= 1 && dT >= 0 && dT <= 1)
1367 {
1368 // Collision detected, calculate intersection coordinates
1369 dXIntersect = dX1 + (dT * dDiffX1);
1370 dYIntersect = dY1 + (dT * dDiffY1);
1371
1372 // And calc the average end-point coordinates
1373 dXAvgEnd = (dX2 + dX4) / 2;
1374 dYAvgEnd = (dY2 + dY4) / 2;
1375
1376 // Get the line segments at which intersection occurred
1377 nProfile1LineSegment = i;
1378 nProfile2LineSegment = j;
1379
1380 // LogStream << "\t" << "INTERSECTION dX2 = " << dX2 << " dX4 = " << dX4 << " dY2 = " << dY2 << " dY4 = " << dY4 << endl;
1381 return true;
1382 }
1383 }
1384 }
1385
1386 // No intersection
1387 return false;
1388}
1389
1390//===============================================================================================================================
1392//===============================================================================================================================
1393void CSimulation::MarkProfilesOnGrid(int const nCoast, int& nValidProfiles)
1394{
1395 // How many profiles on this coast?
1396 int const nProfiles = m_VCoast[nCoast].nGetNumProfiles();
1397
1398 if (nProfiles == 0)
1399 {
1400 // This can happen if the coastline is very short, so just give a warning and carry on with the next coastline
1402 LogStream << WARN << m_ulIter << ": coast " << nCoast << " has no profiles" << endl;
1403
1404 return;
1405 }
1406
1407 static bool bDownCoast = true;
1408
1409 // Now do this for every profile, alternate between up-coast and down-coast directions
1410 for (int n = 0; n < nProfiles; n++)
1411 {
1412 CGeomProfile* pProfile;
1413
1414 if (bDownCoast)
1415 pProfile = m_VCoast[nCoast].pGetProfileWithDownCoastSeq(n);
1416 else
1417 pProfile = m_VCoast[nCoast].pGetProfileWithUpCoastSeq(n);
1418
1419 // Don't do this for the first and last profiles (i.e. the profiles at the start and end of the coast) since these are put onto the grid elsewhere
1420 if ((pProfile->bStartOfCoast()) || (pProfile->bEndOfCoast()))
1421 continue;
1422
1423 int const nProfile = pProfile->nGetProfileID();
1424
1425 // If this profile has a problem, then forget about it
1426 // if (! pProfile->bProfileOK())
1427 // {
1428 // LogStream << m_ulIter << ": in MarkProfilesOnGrid() profile " << nProfile << " is not OK" << endl;
1429 // continue;
1430 // }
1431
1432 int const nPoints = pProfile->nGetProfileSize();
1433
1434 if (nPoints < 2)
1435 {
1436 // Need at least two points in the profile, so this profile is invalid: mark it
1437 m_VCoast[nCoast].pGetProfile(nProfile)->SetTooShort(true);
1438
1440 LogStream << m_ulIter << ": coast " << nCoast << ", profile " << nProfile << " is invalid, has only " << nPoints << " points" << endl;
1441
1442 continue;
1443 }
1444
1445 // OK, go for it: set up temporary vectors to hold the x-y coords (in grid CRS) of the cells which we will mark
1446 vector<CGeom2DIPoint> VCellsToMark;
1447 vector<bool> bVShared;
1448 bool bTooShort = false;
1449 bool bTruncatedSameCoast = false;
1450 bool bHitCoast = false;
1451 bool bHitLand = false;
1452 bool bHitIntervention = false;
1453 bool bHitAnotherProfile = false; // TODO 044
1454
1455 CreateRasterizedProfile(nCoast, pProfile, &VCellsToMark, &bVShared, bTooShort, bTruncatedSameCoast, bHitCoast, bHitLand, bHitIntervention, bHitAnotherProfile); // TODO 044
1456
1457 if ((bTruncatedSameCoast && (! ACCEPT_TRUNCATED_PROFILES)) || bTooShort || bHitCoast || bHitLand || bHitIntervention || bHitAnotherProfile || VCellsToMark.size() == 0)
1458 continue;
1459
1460 // This profile is fine
1461 nValidProfiles++;
1462
1463 for (unsigned int k = 0; k < VCellsToMark.size(); k++)
1464 {
1465 // Ignore duplicate points
1466 if ((k > 0) && (VCellsToMark[k] == m_VCoast[nCoast].pGetProfile(nProfile)->pPtiGetLastCellInProfile()))
1467 continue;
1468
1469 // Mark each cell in the raster grid
1470 int nXTmp = VCellsToMark[k].nGetX();
1471 int nYTmp = VCellsToMark[k].nGetY();
1472 m_pRasterGrid->m_Cell[nXTmp][nYTmp].SetCoastAndProfileID(nCoast, nProfile);
1473
1474 // Store the raster grid coordinates in the profile object
1475 m_VCoast[nCoast].pGetProfile(nProfile)->AppendCellInProfile(nXTmp, nYTmp);
1476
1477 // Mark the shared (i.e. multi-line) parts of the profile (if any)
1478 // if (bVShared[k])
1479 // {
1480 // m_VCoast[nCoast].pGetProfile(nProfile)->AppendPointShared(true);
1481 // // LogStream << m_ulIter << ": profile " << j << " point " << k << " marked as shared" << endl;
1482 // }
1483 // else
1484 // {
1485 // m_VCoast[nCoast].pGetProfile(nProfile)->AppendPointShared(false);
1486 // // LogStream << m_ulIter << ": profile " << nProfile << " point " << k << " marked as NOT shared" << endl;
1487 // }
1488 }
1489
1490 // Get the deep water wave height and orientation values at the end of the profile
1491 double const dDeepWaterWaveHeight = m_pRasterGrid->m_Cell[VCellsToMark.back().nGetX()][VCellsToMark.back().nGetY()].dGetCellDeepWaterWaveHeight();
1492 double const dDeepWaterWaveAngle = m_pRasterGrid->m_Cell[VCellsToMark.back().nGetX()][VCellsToMark.back().nGetY()].dGetCellDeepWaterWaveAngle();
1493 double const dDeepWaterWavePeriod = m_pRasterGrid->m_Cell[VCellsToMark.back().nGetX()][VCellsToMark.back().nGetY()].dGetCellDeepWaterWavePeriod();
1494
1495 // And store them for this profile
1496 m_VCoast[nCoast].pGetProfile(nProfile)->SetProfileDeepWaterWaveHeight(dDeepWaterWaveHeight);
1497 m_VCoast[nCoast].pGetProfile(nProfile)->SetProfileDeepWaterWaveAngle(dDeepWaterWaveAngle);
1498 m_VCoast[nCoast].pGetProfile(nProfile)->SetProfileDeepWaterWavePeriod(dDeepWaterWavePeriod);
1499 }
1500
1501 bDownCoast = ! bDownCoast;
1502}
1503
1504//===============================================================================================================================
1506//===============================================================================================================================
1507void CSimulation::CreateRasterizedProfile(int const nCoast, CGeomProfile* pProfile, vector<CGeom2DIPoint>* pVIPointsOut, vector<bool>* pbVShared, bool& bTooShort, bool& bTruncatedSameCoast, bool& bHitCoast, bool& bHitLand, bool& bHitIntervention, bool& bHitAnotherProfile) // TODO 044
1508{
1509 int const nProfile = pProfile->nGetProfileID();
1510 int nSeg = 0;
1511 int const nNumSegments = pProfile->nGetNumLineSegments();
1512
1513 pVIPointsOut->clear();
1514
1515 // LogStream << m_ulIter << ": in CreateRasterizedProfile() *pPtiStart for profile " << nProfile << " is [" << pPtiStart->nGetX() << "][" << pPtiStart->nGetY() << "]" << endl;
1516 int nXStartLast = INT_NODATA;
1517 int nYStartLast = INT_NODATA;
1518 int nXEndLast = INT_NODATA;
1519 int nYEndLast = INT_NODATA;
1520
1521 // Do for every segment of this profile
1522 for (nSeg = 0; nSeg < nNumSegments; nSeg++)
1523 {
1524 // Do once for every line segment
1525 CGeom2DIPoint PtiSegStart;
1526
1527 if (nSeg == 0)
1528 {
1529 // If this is the first segment, use the profile start point to prevent Ext CRS to Grid CGS rounding errors
1530 PtiSegStart = pProfile->pPtiGetStartPoint();
1531 }
1532 else
1533 {
1534 CGeom2DPoint const* pPtSegStart = pProfile->pPtGetPointInProfile(nSeg);
1535
1536 // Convert from the external CRS to grid CRS
1537 PtiSegStart = PtiExtCRSToGridRound(pPtSegStart);
1538 }
1539
1540 CGeom2DPoint const* pPtSegEnd = pProfile->pPtGetPointInProfile(nSeg + 1); // This is OK
1541
1542 // Convert from the external CRS to grid CRS
1543 CGeom2DIPoint const PtiSegEnd = PtiExtCRSToGridRound(pPtSegEnd);
1544
1545 // Safety check
1546 if (PtiSegStart == PtiSegEnd)
1547 continue;
1548
1549 int const nXStart = PtiSegStart.nGetX();
1550 int const nYStart = PtiSegStart.nGetY();
1551 int const nXEnd = PtiSegEnd.nGetX();
1552 int const nYEnd = PtiSegEnd.nGetY();
1553
1554 bool bShared = false;
1555
1556 if (pProfile->nGetNumCoincidentProfilesInLineSegment(nSeg) > 1)
1557 {
1558 bShared = true;
1559
1560 // If this is the second or more of several coincident line segments (i.e. it has the same start and end points as the previous line segment) then ignore it
1561 if ((nXStart == nXStartLast) && (nYStart == nYStartLast) && (nXEnd == nXEndLast) && (nYEnd == nYEndLast))
1562 continue;
1563 }
1564
1565 // Interpolate between cells by a simple DDA line algorithm, see http://en.wikipedia.org/wiki/Digital_differential_analyzer_(graphics_algorithm) Note that Bresenham's algorithm gave occasional gaps
1566 double dXInc = nXEnd - nXStart;
1567 double dYInc = nYEnd - nYStart;
1568 double const dLength = tMax(tAbs(dXInc), tAbs(dYInc));
1569
1570 dXInc /= dLength;
1571 dYInc /= dLength;
1572
1573 double dX = nXStart;
1574 double dY = nYStart;
1575
1576 // Process each interpolated point
1577 for (int m = 0; m <= nRound(dLength); m++)
1578 {
1579 int const nX = nRound(dX);
1580 int const nY = nRound(dY);
1581
1582 // Do some checking of this interpolated point, but only if this is not a grid-edge profile (these profiles are always valid)
1583 if ((! pProfile->bStartOfCoast()) && (! pProfile->bEndOfCoast()))
1584 {
1585 // Is the interpolated point within the valid raster grid?
1586 if (! bIsWithinValidGrid(nX, nY))
1587 {
1588 // It is outside the valid grid, so mark this profile and quit the loop
1589 bTruncatedSameCoast = true;
1590
1592 pProfile->SetTruncatedSameCoast(true);
1593
1595 LogStream << m_ulIter << ": profile " << nProfile << " is invalid, truncated at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
1596
1597 break;
1598 }
1599
1600 // Check again: is this cell (or an adjacent cell: does not matter which) already marked as 'under' a profile? TODO check that nY+1 is OK
1601 if (m_pRasterGrid->m_Cell[nX][nY].bIsProfile() || m_pRasterGrid->m_Cell[nX][nY+1].bIsProfile())
1602 {
1603 // This cell or an adjacent cell, is 'under' a profile, so now check if the profile belongs to another coast
1604 int nHitProfileCoast1 = m_pRasterGrid->m_Cell[nX][nY].nGetProfileCoastID();
1605 int nHitProfileCoast2 = m_pRasterGrid->m_Cell[nX][nY+1].nGetProfileCoastID();
1606
1607 if ((nHitProfileCoast1 == nCoast) || (nHitProfileCoast2 == nCoast))
1608 {
1609 // The profile belongs to the same coast, mark this profile as invalid
1610 bHitAnotherProfile = true;
1611 pProfile->SetHitAnotherProfile(true);
1612 return;
1613 }
1614 }
1615
1616 // If this is the first line segment of the profile, then once we are clear of the coastline (say, when m > 3), check if this profile hits land at this interpolated point. NOTE Get problems here since if the coastline vector has been heavily smoothed, this can result is 'false positives' profiles marked as invalid which are not actually invalid, because the profile hits land when m = 0 or m = 1. This results in some cells being flagged as profile cells which are actually inland
1618 {
1619 // Check this cell. Two diagonal(ish) raster lines can cross each other without any intersection, so must also test an adjacent cell for intersection (does not matter which adjacent cell)
1620 if ((m_pRasterGrid->m_Cell[nX][nY].bIsCoastline()) || (bIsWithinValidGrid(nX, nY + 1) && m_pRasterGrid->m_Cell[nX][nY + 1].bIsCoastline()))
1621 {
1622 // We've hit a coastline so set a switch and mark the profile, then quit
1623 bHitCoast = true;
1624 pProfile->SetHitCoast(true);
1625 int const nHitCoast = m_pRasterGrid->m_Cell[nX][nY].nGetCoastline();
1626
1628 LogStream << m_ulIter << ": coast " << nCoast << " profile " << nProfile << " is invalid, hit coast " << nHitCoast << " at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
1629
1630 return;
1631 }
1632
1633 if (! m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea())
1634 {
1635 // We've hit dry land, so set a switch and mark the profile
1636 bHitLand = true;
1637 pProfile->SetHitLand(true);
1638
1639 LogStream << m_ulIter << ": profile " << nProfile << " HIT LAND at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}, elevation = " << m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() << ", SWL = " << m_dThisIterSWL << endl;
1640
1641 return;
1642 }
1643
1644 if (m_pRasterGrid->m_Cell[nX][nY].nGetInterventionClass() != INT_NODATA)
1645 {
1646 // We've hit an intervention, so set a switch and mark the profile
1647 bHitIntervention = true;
1648 pProfile->SetHitIntervention(true);
1649
1650 LogStream << m_ulIter << ": profile " << nProfile << " HIT INTERVENTION at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}, elevation = " << m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() << ", SWL = " << m_dThisIterSWL << endl;
1651
1652 return;
1653 }
1654 }
1655
1656 // Now check to see if we hit another profile which is not a coincident normal to this normal
1657 if (m_pRasterGrid->m_Cell[nX][nY].bIsProfile())
1658 {
1659 // We've hit a raster cell which is already marked as 'under' a normal profile. Get the number of the profile which marked this cell, and the coast to hich this profile belongs
1660 int const nHitProfile = m_pRasterGrid->m_Cell[nX][nY].nGetProfileID();
1661 int const nHitProfileCoast = m_pRasterGrid->m_Cell[nX][nY].nGetProfileCoastID();
1662
1663 // Do both profiles belong to the same coast?
1664 if (nCoast == nHitProfileCoast)
1665 {
1666 // Both profiles belong to the same coast. Is this the number of a coincident profile of this profile?
1667 if (! pProfile->bFindProfileInCoincidentProfilesOfLastLineSegment(nHitProfile))
1668 {
1669 // It isn't a coincident profile, so we have just hit an unrelated profile. Mark this profile as invalid and move on
1670 pProfile->SetHitAnotherProfile(true);
1671 bHitAnotherProfile = true;
1672
1673 return;
1674 }
1675 }
1676 }
1677 }
1678
1679 // Append this point to the output vector
1680 pVIPointsOut->push_back(CGeom2DIPoint(nX, nY)); // Is in raster grid coordinates
1681 pbVShared->push_back(bShared);
1682
1683 // And increment for next time
1684 dX += dXInc;
1685 dY += dYInc;
1686 }
1687
1688 nXStartLast = nXStart;
1689 nYStartLast = nYStart;
1690 nXEndLast = nXEnd;
1691 nYEndLast = nYEnd;
1692
1693 if (bTruncatedSameCoast)
1694 break;
1695 }
1696
1697 if (bTruncatedSameCoast)
1698 {
1699 if (nSeg < (nNumSegments - 1))
1700 // We are truncating the profile, so remove any line segments after this one
1701 pProfile->TruncateLineSegments(nSeg);
1702
1703 // Shorten the vector input. Ignore CPPCheck errors here, since we know that pVIPointsOut is not empty
1704 int const nLastX = pVIPointsOut->at(pVIPointsOut->size() - 1).nGetX();
1705 int const nLastY = pVIPointsOut->at(pVIPointsOut->size() - 1).nGetY();
1706
1707 pProfile->pPtGetPointInProfile(nSeg + 1)->SetX(dGridCentroidXToExtCRSX(nLastX));
1708 pProfile->pPtGetPointInProfile(nSeg + 1)->SetY(dGridCentroidYToExtCRSY(nLastY));
1709 }
1710
1711 // // DEBUG CODE =====================================================================================
1712 // LogStream << "====================" << endl;
1713 // LogStream << m_ulIter << ": for profile " << nProfile << " pPtiStart = [" << pPtiStart->nGetX() << "][" << pPtiStart->nGetY() << "] pPtiEnd = [" << pPtiEnd->nGetX() << "][" << pPtiEnd->nGetY() << "] pVIPointsOut->size() = " << pVIPointsOut->size() << endl;
1714 // // for (int n = 0; n < static_cast<int>(pVIPointsOut->size()); n++)
1715 // // LogStream << "\t[" << pVIPointsOut->at(n).nGetX() << "][" << pVIPointsOut->at(n).nGetY() << "]" << endl;
1716 // LogStream << "====================" << endl;
1717 // // DEBUG CODE =====================================================================================
1718
1719 if (pVIPointsOut->size() < 3)
1720 {
1721 // Coastline-normal profiles cannot be very short (e.g. with less than 3 cells), since we cannot calculate along-profile slope properly for such short profiles
1722 bTooShort = true;
1723 pProfile->SetTooShort(true);
1724
1726 {
1727 // Ignore CPPCheck errors here, since we know that pVIPointsOut is not empty
1728 LogStream << m_ulIter << ": profile " << nProfile << " is invalid, is too short, only " << pVIPointsOut->size() << " points, HitLand?" << bHitLand << ". From [" << pVIPointsOut->at(0).nGetX() << "][" << pVIPointsOut->at(0).nGetY() << "] = {" << dGridCentroidXToExtCRSX(pVIPointsOut->at(0).nGetX()) << ", " << dGridCentroidYToExtCRSY(pVIPointsOut->at(0).nGetY()) << "} to [" << pVIPointsOut->at(pVIPointsOut->size() - 1).nGetX() << "][" << pVIPointsOut->at(pVIPointsOut->size() - 1).nGetY() << "] = {" << dGridCentroidXToExtCRSX(pVIPointsOut->at(pVIPointsOut->size() - 1).nGetX()) << ", " << dGridCentroidYToExtCRSY(pVIPointsOut->at(pVIPointsOut->size() - 1).nGetY()) << "}" << endl;
1729 }
1730 }
1731}
1732
1733//===============================================================================================================================
1735//===============================================================================================================================
1736void CSimulation::MergeProfilesAtFinalLineSegments(int const nCoast, CGeomProfile* pFirstProfile, CGeomProfile* pSecondProfile, int const nFirstProfileLineSegments, int const nSecondProfileLineSegments, double const dIntersectX, double const dIntersectY, double const dAvgEndX, double const dAvgEndY)
1737{
1738 // The point of intersection is on the final (most seaward) line segment of both profiles. Put together a vector of coincident profile numbers (with no duplicates) for both profiles
1739 int nCombinedLastSeg = 0;
1740 vector<pair<int, int>> prVCombinedProfilesCoincidentProfilesLastSeg;
1741
1742 for (unsigned int n = 0; n < pFirstProfile->pprVGetPairedCoincidentProfilesForLineSegment(nFirstProfileLineSegments - 1)->size(); n++)
1743 {
1744 pair<int, int> prTmp;
1745 prTmp.first = pFirstProfile->pprVGetPairedCoincidentProfilesForLineSegment(nFirstProfileLineSegments - 1)->at(n).first;
1746 prTmp.second = pFirstProfile->pprVGetPairedCoincidentProfilesForLineSegment(nFirstProfileLineSegments - 1)->at(n).second;
1747
1748 bool bFound = false;
1749
1750 for (unsigned int m = 0; m < prVCombinedProfilesCoincidentProfilesLastSeg.size(); m++)
1751 {
1752 if (prVCombinedProfilesCoincidentProfilesLastSeg[m].first == prTmp.first)
1753 {
1754 bFound = true;
1755 break;
1756 }
1757 }
1758
1759 if (! bFound)
1760 {
1761 prVCombinedProfilesCoincidentProfilesLastSeg.push_back(prTmp);
1762 nCombinedLastSeg++;
1763 }
1764 }
1765
1766 for (unsigned int n = 0; n < pSecondProfile->pprVGetPairedCoincidentProfilesForLineSegment(nSecondProfileLineSegments - 1)->size(); n++)
1767 {
1768 pair<int, int> prTmp;
1769 prTmp.first = pSecondProfile->pprVGetPairedCoincidentProfilesForLineSegment(nSecondProfileLineSegments - 1)->at(n).first;
1770 prTmp.second = pSecondProfile->pprVGetPairedCoincidentProfilesForLineSegment(nSecondProfileLineSegments - 1)->at(n).second;
1771
1772 bool bFound = false;
1773
1774 for (unsigned int m = 0; m < prVCombinedProfilesCoincidentProfilesLastSeg.size(); m++)
1775 {
1776 if (prVCombinedProfilesCoincidentProfilesLastSeg[m].first == prTmp.first)
1777 {
1778 bFound = true;
1779 break;
1780 }
1781 }
1782
1783 if (! bFound)
1784 {
1785 prVCombinedProfilesCoincidentProfilesLastSeg.push_back(prTmp);
1786 nCombinedLastSeg++;
1787 }
1788 }
1789
1790 // Increment the number of each line segment
1791 for (int m = 0; m < nCombinedLastSeg; m++)
1792 prVCombinedProfilesCoincidentProfilesLastSeg[m].second++;
1793
1794 vector<pair<int, int>> prVFirstProfileCoincidentProfilesLastSeg = *pFirstProfile->pprVGetPairedCoincidentProfilesForLineSegment(nFirstProfileLineSegments - 1);
1795 vector<pair<int, int>> prVSecondProfileCoincidentProfilesLastSeg = *pSecondProfile->pprVGetPairedCoincidentProfilesForLineSegment(nSecondProfileLineSegments - 1);
1796 int const nNumFirstProfileCoincidentProfilesLastSeg = static_cast<int>(prVFirstProfileCoincidentProfilesLastSeg.size());
1797 int const nNumSecondProfileCoincidentProfilesLastSeg = static_cast<int>(prVSecondProfileCoincidentProfilesLastSeg.size());
1798
1799 // LogStream << m_ulIter << ": END-SEGMENT INTERSECTION between profiles " << nFirstProfile << " and " << nSecondProfile << " at line segment " << nFirstProfileLineSegments-1 << "/" << nFirstProfileLineSegments-1 << ", and line segment " << nSecondProfileLineSegments-1 << "/" << nSecondProfileLineSegments-1 << ", respectively. Both truncated at [" << dIntersectX << ", " << dIntersectY << "] then profiles {" << nFirstProfile << "} and {" << nSecondProfile << "} extended to [" << dAvgEndX << ", " << dAvgEndY << "]" << endl;
1800
1801 // Truncate the first profile, and all co-incident profiles, at the point of intersection
1802 for (int n = 0; n < nNumFirstProfileCoincidentProfilesLastSeg; n++)
1803 {
1804 int const nThisProfile = prVFirstProfileCoincidentProfilesLastSeg[n].first;
1805 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
1806 int const nProfileLength = pThisProfile->nGetProfileSize();
1807
1808 // This is the final line segment of the first 'main' profile. We are assuming that it is also the final line segment of all co-incident profiles. This is fine, altho' each profile may well have a different number of line segments landwards i.e. the number of the line segment may be different for each co-incident profile
1809 pThisProfile->SetPointInProfile(nProfileLength - 1, dIntersectX, dIntersectY);
1810 }
1811
1812 // Truncate the second profile, and all co-incident profiles, at the point of intersection
1813 for (int n = 0; n < nNumSecondProfileCoincidentProfilesLastSeg; n++)
1814 {
1815 int const nThisProfile = prVSecondProfileCoincidentProfilesLastSeg[n].first;
1816 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
1817 int const nProfileLength = pThisProfile->nGetProfileSize();
1818
1819 // This is the final line segment of the second 'main' profile. We are assuming that it is also the final line segment of all co-incident profiles. This is fine, altho' each profile may well have a different number of line segments landwards i.e. the number of the line segment may be different for each co-incident profile
1820 pThisProfile->SetPointInProfile(nProfileLength - 1, dIntersectX, dIntersectY);
1821 }
1822
1823 // Append a new straight line segment to the existing line segment(s) of the first profile, and to all co-incident profiles
1824 for (int nThisLineSeg = 0; nThisLineSeg < nNumFirstProfileCoincidentProfilesLastSeg; nThisLineSeg++)
1825 {
1826 int const nThisProfile = prVFirstProfileCoincidentProfilesLastSeg[nThisLineSeg].first;
1827 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
1828
1829 // Update this profile
1830 pThisProfile->AppendPointInProfile(dAvgEndX, dAvgEndY);
1831
1832 // Append details of the combined profiles
1833 pThisProfile->AppendLineSegment();
1834
1835 for (int m = 0; m < nCombinedLastSeg; m++)
1836 pThisProfile->AppendCoincidentProfileToLineSegments(prVCombinedProfilesCoincidentProfilesLastSeg[m]);
1837 }
1838
1839 // Append a new straight line segment to the existing line segment(s) of the second profile, and to all co-incident profiles
1840 for (int nThisLineSeg = 0; nThisLineSeg < nNumSecondProfileCoincidentProfilesLastSeg; nThisLineSeg++)
1841 {
1842 int const nThisProfile = prVSecondProfileCoincidentProfilesLastSeg[nThisLineSeg].first;
1843 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
1844
1845 // Update this profile
1846 pThisProfile->AppendPointInProfile(dAvgEndX, dAvgEndY);
1847
1848 // Append details of the combined profiles
1849 pThisProfile->AppendLineSegment();
1850
1851 for (int m = 0; m < nCombinedLastSeg; m++)
1852 pThisProfile->AppendCoincidentProfileToLineSegments(prVCombinedProfilesCoincidentProfilesLastSeg[m]);
1853 }
1854
1855 // // DEBUG CODE ****************************************************************
1856 // int nFirstProfileLineSeg= pFirstProfile->nGetNumLineSegments();
1857 // int nSecondProfileLineSeg = pSecondProfile->nGetNumLineSegments();
1858 //
1859 // LogStream << "\tProfile {" << nFirstProfile << "} now has " << nFirstProfileLineSeg << " line segments" << endl;
1860 // for (int m = 0; m < nFirstProfileLineSeg; m++)
1861 // {
1862 // vector<pair<int, int> > prVCoincidentProfiles = *pFirstProfile->pprVGetPairedCoincidentProfilesForLineSegment(m);
1863 // LogStream << "\tCo-incident profiles and line segments for line segment " << m << " of profile {" << nFirstProfile << "} are {";
1864 // for (int nn = 0; nn < prVCoincidentProfiles.size(); nn++)
1865 // LogStream << " " << prVCoincidentProfiles[nn].first << "[" << prVCoincidentProfiles[nn].second << "] ";
1866 // LogStream << " }" << endl;
1867 // }
1868 // LogStream << "\tProfile {" << nSecondProfile << "} now has " << nSecondProfileLineSeg << " line segments" << endl;
1869 // for (int m = 0; m < nSecondProfileLineSeg; m++)
1870 // {
1871 // vector<pair<int, int> > prVCoincidentProfiles = *pSecondProfile->pprVGetPairedCoincidentProfilesForLineSegment(m);
1872 // LogStream << "\tCo-incident profiles and line segments for line segment " << m << " of profile {" << nSecondProfile << "} are {";
1873 // for (int nn = 0; nn < prVCoincidentProfiles.size(); nn++)
1874 // LogStream << " " << prVCoincidentProfiles[nn].first << "[" << prVCoincidentProfiles[nn].second << "] ";
1875 // LogStream << " }" << endl;
1876 // }
1877 // // DEBUG CODE ******************************************************************
1878}
1879
1880//===============================================================================================================================
1882//===============================================================================================================================
1883void CSimulation::TruncateOneProfileRetainOtherProfile(int const nCoast, CGeomProfile* pProfileToTruncate, CGeomProfile* pProfileToRetain, double dIntersectX, double dIntersectY, int nProfileToTruncateIntersectLineSeg, int nProfileToRetainIntersectLineSeg, bool const bAlreadyPresent)
1884{
1885 // // Occasionally, profiles cross each other, with the crossing not detected. So check for intersection between pProfileToTruncate and all profiles (starting from the last i.e. in reverse order) in all segments (starting from the last i.e. in reverse order) of pProfileToRetain's CGeomMultiLine
1886 // bool bFound = false;
1887 // int nNumSegProfRetain = pProfileToRetain->nGetNumLineSegments();
1888 // for (int nSeg = nNumSegProfRetain-1; nSeg >= 0; nSeg--)
1889 // {
1890 // if (bFound)
1891 // break;
1892 //
1893 // int nNumProfInSeg = pProfileToRetain->nGetNumCoincidentProfilesInLineSegment(nSeg);
1894 // for (int nProf = nNumProfInSeg-1; nProf >= 0; nProf--)
1895 // {
1896 // int nThisProf = pProfileToRetain->nGetCoincidentProfileForLineSegment(nSeg, nProf);
1897 // CGeomProfile* pThisProf = m_VCoast[nCoast].pGetProfile(nThisProf);
1898 //
1899 // int nProfToTruncLineSeg = 0;
1900 // int nThisProfLineSeg = 0;
1901 // double dTmpIntersectX = 0;
1902 // double dTmpIntersectY = 0;
1903 // double dAvgEndX = 0;
1904 // double dAvgEndY = 0;
1905 //
1906 // if (bCheckForIntersection(pProfileToTruncate, pThisProf, nProfToTruncLineSeg, nThisProfLineSeg, dTmpIntersectX, dTmpIntersectY, dAvgEndX, dAvgEndY))
1907 // {
1908 // // An intersection was found: so the profile with which pProfileToTruncate intersects becomes the new pProfileToRetain, and dIntersectX, dIntersectY, nProfileToTruncateIntersectLineSeg, and nProfileToRetainIntersectLineSeg are also changed
1909 // pProfileToRetain = pThisProf;
1910 // dIntersectX = dTmpIntersectX;
1911 // dIntersectY = dTmpIntersectY;
1912 // nProfileToRetainIntersectLineSeg = nThisProfLineSeg;
1913 // nProfileToTruncateIntersectLineSeg = nProfToTruncLineSeg;
1914 //
1915 // bFound = true;
1916 // break;
1917 // }
1918 // }
1919 // }
1920
1921 // Insert the intersection point into the main retain-profile if it is not already in the profile, and do the same for all co-incident profiles of the main retain-profile. Also add details of the to-truncate profile (and all its coincident profiles) to every line segment of the main to-retain profile which is seaward of the point of intersection
1922 int const nRet = nInsertPointIntoProfilesIfNeededThenUpdate(nCoast, pProfileToRetain, dIntersectX, dIntersectY, nProfileToRetainIntersectLineSeg, pProfileToTruncate, nProfileToTruncateIntersectLineSeg, bAlreadyPresent);
1923
1924 if (nRet != RTN_OK)
1925 {
1926 // LogStream << m_ulIter << ": error in nInsertPointIntoProfilesIfNeededThenUpdate()" << endl;
1927 return;
1928 }
1929
1930 // Get all profile points of the main retain-profile seawards from the intersection point, and do the same for the corresponding line segments (including coincident profiles). This also includes details of the main to-truncate profile (and all its coincident profiles)
1931 vector<CGeom2DPoint> PtVProfileLastPart;
1932 vector<vector<pair<int, int>>> prVLineSegLastPart;
1933
1934 if (bAlreadyPresent)
1935 {
1936 PtVProfileLastPart = pProfileToRetain->PtVGetThisPointAndAllAfter(nProfileToRetainIntersectLineSeg);
1937 prVLineSegLastPart = pProfileToRetain->prVVGetAllLineSegAfter(nProfileToRetainIntersectLineSeg);
1938 }
1939
1940 else
1941 {
1942 PtVProfileLastPart = pProfileToRetain->PtVGetThisPointAndAllAfter(nProfileToRetainIntersectLineSeg + 1);
1943 prVLineSegLastPart = pProfileToRetain->prVVGetAllLineSegAfter(nProfileToRetainIntersectLineSeg + 1);
1944 }
1945
1946 // assert(PtVProfileLastPart.size() > 1);
1947 // assert(prVLineSegLastPart.size() > 0);
1948
1949 // Truncate the truncate-profile at the point of intersection, and do the same for all its co-incident profiles. Then append the profile points of the main to-retain profile seaward from the intersection point, and do the same for the corresponding line segments (including coincident profiles)
1950 TruncateProfileAndAppendNew(nCoast, pProfileToTruncate, nProfileToTruncateIntersectLineSeg, &PtVProfileLastPart, &prVLineSegLastPart);
1951
1952 // assert(m_VCoast[nCoast].pGetProfile(nProfileToTruncate)->nGetProfileSize() > 1);
1953 // assert(pProfileToRetain->nGetNumLineSegments() > 0);
1954 // assert(m_VCoast[nCoast].pGetProfile(nProfileToTruncate)->nGetNumLineSegments() > 0);
1955}
1956
1957//===============================================================================================================================
1959//===============================================================================================================================
1960int CSimulation::nInsertPointIntoProfilesIfNeededThenUpdate(int const nCoast, CGeomProfile* pProfileToRetain, double const dIntersectX, double const dIntersectY, int const nProfileToRetainIntersectLineSeg, CGeomProfile* pProfileToTruncate, int const nProfileToTruncateIntersectLineSeg, bool const bAlreadyPresent)
1961{
1962 // // DEBUG CODE ****************************************************************
1963 // // Get the index numbers of all coincident profiles for the 'main' to-retain profile for the line segment in which intersection occurred
1964 // vector<pair<int, int> > prVRetainCoincidentProfilesCHECK1 = *m_VCoast[nCoast].pGetProfile(nMainProfile)->pprVGetPairedCoincidentProfilesForLineSegment(nProfileToRetainIntersectLineSeg);
1965 // int nNumRetainCoincidentCHECK1 = prVRetainCoincidentProfilesCHECK1.size();
1966 // for (int nn = 0; nn < nNumRetainCoincidentCHECK1; nn++)
1967 // {
1968 // int nThisProfile = prVRetainCoincidentProfilesCHECK1[nn].first;
1969 // LogStream << "\tBEFORE nInsertPointIntoProfilesIfNeededThenUpdate(): " << (nThisProfile == nMainProfile ? "MAIN" : "COINCIDENT") << " to-retain profile {" << nThisProfile << "} has " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize() << " points ";
1970 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize(); nn++)
1971 // LogStream << "[" << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetX() << ", " << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetY() << "] ";
1972 // LogStream << "), and " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments() << " line segments, co-incident profiles and their line segments are ";
1973 // for (int mm = 0; mm < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments(); mm++)
1974 // {
1975 // vector<pair<int, int> > prVTmp = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pprVGetPairedCoincidentProfilesForLineSegment(mm);
1976 // LogStream << "{ ";
1977 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumCoincidentProfilesInLineSegment(mm); nn++)
1978 // LogStream << prVTmp[nn].first << "[" << prVTmp[nn].second << "] ";
1979 // LogStream << "} ";
1980 // }
1981 // LogStream << endl;
1982 //
1983 // for (int nPoint = 0; nPoint < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize()-1; nPoint++)
1984 // {
1985 // CGeom2DPoint
1986 // Pt1 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint),
1987 // Pt2 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint+1);
1988 //
1989 // if (Pt1 == Pt2)
1990 // LogStream << m_ulIter << ": IDENTICAL POINTS before changes, in profile {" << nThisProfile << "} points = " << nPoint << " and " << nPoint+1 << endl;
1991 // }
1992 // }
1993 //
1994 // // Get the index numbers of all coincident profiles for the 'main' to-truncate profile for the line segment in which intersection occurred
1995 // vector<pair<int, int> > prVTruncateCoincidentProfilesCHECK1 = *m_VCoast[nCoast].pGetProfile(nProfileToTruncate)->pprVGetPairedCoincidentProfilesForLineSegment(nProfileToTruncateIntersectLineSeg);
1996 // int nNumTruncateCoincidentCHECK1 = prVTruncateCoincidentProfilesCHECK1.size();
1997 // for (int nn = 0; nn < nNumTruncateCoincidentCHECK1; nn++)
1998 // {
1999 // int nThisProfile = prVTruncateCoincidentProfilesCHECK1[nn].first;
2000 // LogStream << "\tBEFORE nInsertPointIntoProfilesIfNeededThenUpdate(): " << (nThisProfile == nProfileToTruncate ? "MAIN" : "COINCIDENT") << " to-truncate profile {" << nThisProfile << "} has " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize() << " points ";
2001 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize(); nn++)
2002 // LogStream << "[" << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetX() << ", " << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetY() << "] ";
2003 // LogStream << "), and " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments() << " line segments, co-incident profiles and their line segments are ";
2004 // for (int mm = 0; mm < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments(); mm++)
2005 // {
2006 // vector<pair<int, int> > prVTmp = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pprVGetPairedCoincidentProfilesForLineSegment(mm);
2007 // LogStream << "{ ";
2008 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumCoincidentProfilesInLineSegment(mm); nn++)
2009 // LogStream << prVTmp[nn].first << "[" << prVTmp[nn].second << "] ";
2010 // LogStream << "} ";
2011 // }
2012 // LogStream << endl;
2013 //
2014 // for (int nPoint = 0; nPoint < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize()-1; nPoint++)
2015 // {
2016 // CGeom2DPoint
2017 // Pt1 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint),
2018 // Pt2 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint+1);
2019 //
2020 // if (Pt1 == Pt2)
2021 // LogStream << m_ulIter << ": IDENTICAL POINTS before changes, in profile {" << nThisProfile << "} points = " << nPoint << " and " << nPoint+1 << endl;
2022 // }
2023 // }
2024 // // DEBUG CODE ******************************************************************
2025
2026 int const nProfileToRetain = pProfileToRetain->nGetProfileID();
2027
2028 // Get the index numbers of all coincident profiles for the 'main' to-retain profile for the line segment in which intersection occurs
2029 vector<pair<int, int>> prVCoincidentProfiles = *pProfileToRetain->pprVGetPairedCoincidentProfilesForLineSegment(nProfileToRetainIntersectLineSeg);
2030 int const nNumCoincident = static_cast<int>(prVCoincidentProfiles.size());
2031 vector<int> nLineSegAfterIntersect(nNumCoincident, -1); // The line segment after the point of intersection, for each co-incident profile
2032
2033 // Do this for the main profile and all profiles which are co-incident for this line segment
2034 for (int nn = 0; nn < nNumCoincident; nn++)
2035 {
2036 int const nThisProfile = prVCoincidentProfiles[nn].first; // The number of this profile
2037 int const nThisLineSeg = prVCoincidentProfiles[nn].second; // The line segment of this profile
2038 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
2039
2040 // Is the intersection point already present in the to-retain profile?
2041 if (! bAlreadyPresent)
2042 {
2043 // It is not already present, so insert it and also update the associated multi-line
2044 if (! pThisProfile->bInsertIntersection(dIntersectX, dIntersectY, nThisLineSeg))
2045 {
2046 // Error
2047 LogStream << WARN << m_ulIter << ": cannot insert a line segment after the final line segment (" << nThisLineSeg << ") for " << (nThisProfile == nProfileToRetain ? "main" : "co-incident") << " profile (" << nThisProfile << "), abandoning" << endl;
2048
2050 }
2051
2052 // LogStream << "\tIntersection point NOT already in " << (nThisProfile == nProfileToRetain ? "main" : "co-incident") << " profile {" << nThisProfile << "}, inserted it as point " << nThisLineSeg+1 << endl;
2053 }
2054
2055 // Get the line segment after intersection
2056 nLineSegAfterIntersect[nn] = nThisLineSeg + 1;
2057 }
2058
2059 // for (int nn = 0; nn < nNumCoincident; nn++)
2060 // LogStream << "\tFor profile " << prVCoincidentProfiles[nn].first << " line segment [" << nLineSegAfterIntersect[nn] << "] is immediately after the intersection point" << endl;
2061
2062 // Get the coincident profiles for the to-truncate profile, at the line segment where intersection occurs
2063 vector<pair<int, int>> prVToTruncateCoincidentProfiles = *pProfileToTruncate->pprVGetPairedCoincidentProfilesForLineSegment(nProfileToTruncateIntersectLineSeg);
2064 int const nNumToTruncateCoincident = static_cast<int>(prVToTruncateCoincidentProfiles.size());
2065
2066 // Now add the number of the to-truncate profile, and all its coincident profiles, to all line segments which are seaward of the point of intersection. Do this for the main profile and all profiles which are co-incident for this line segment
2067 for (int nn = 0; nn < nNumCoincident; nn++)
2068 {
2069 int const nThisProfile = prVCoincidentProfiles[nn].first; // The number of this profile
2070 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
2071
2072 // Get the number of line segments for this to-retain profile (will have just increased, if we just inserted a point)
2073 int const nNumLineSegs = pThisProfile->nGetNumLineSegments();
2074
2075 // Do for all line segments seaward of the point of intersection
2076 for (int nLineSeg = nLineSegAfterIntersect[nn], nIncr = 0; nLineSeg < nNumLineSegs; nLineSeg++, nIncr++)
2077 {
2078 // // This can happen occasionally
2079 // if (nThisProfile == nProfileToTruncateIntersectLineSeg)
2080 // {
2081 // LogStream << "\t*** ERROR nThisProfile = " << nThisProfile << " nProfileToTruncateIntersectLineSeg = " << nProfileToTruncateIntersectLineSeg << ", ignoring" << endl;
2082 // pThisProfile->SetHitAnotherProfile(true);
2083 // continue;
2084 // }
2085
2086 // Add the number of the to-truncate profile, and all its coincident profiles, to this line segment
2087 for (int m = 0; m < nNumToTruncateCoincident; m++)
2088 {
2089 int const nProfileToAdd = prVToTruncateCoincidentProfiles[m].first;
2090 int const nProfileToAddLineSeg = prVToTruncateCoincidentProfiles[m].second;
2091
2092 // LogStream << "\tAdding " << (nProfileToAdd == nProfileToTruncateIntersectLineSeg ? "main" : "co-incident") << " truncate-profile " << nProfileToAdd << ", line segment [" << nProfileToAddLineSeg + nIncr << "] to line segment " << nLineSeg << " of " << (nThisProfile == nProfileToRetain ? "main" : "co-incident") << " to-retain profile " << nThisProfile << endl;
2093
2094 pThisProfile->AddCoincidentProfileToExistingLineSegment(nLineSeg, nProfileToAdd, nProfileToAddLineSeg + nIncr);
2095 }
2096 }
2097 }
2098
2099 // // DEBUG CODE ****************************************************************
2100 // Get the index numbers of all coincident profiles for the 'main' profile for the line segment in which intersection occurred
2101 // vector<pair<int, int> > prVCoincidentProfilesCHECK2 = *m_VCoast[nCoast].pGetProfile(nProfileToRetain)->pprVGetPairedCoincidentProfilesForLineSegment(nProfileToRetainIntersectLineSeg);
2102 // int nNumCoincidentCHECK2 = prVCoincidentProfilesCHECK2.size();
2103 // for (int nn = 0; nn < nNumCoincidentCHECK2; nn++)
2104 // {
2105 // int nThisProfile = prVCoincidentProfilesCHECK2[nn].first;
2106 // LogStream << "\tAFTER nInsertPointIntoProfilesIfNeededThenUpdate(): " << (nThisProfile == nProfileToRetain ? "MAIN" : "COINCIDENT") << " to-retain profile {" << nThisProfile << "} has " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize() << " points ";
2107 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize(); nn++)
2108 // LogStream << "[" << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetX() << ", " << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetY() << "] ";
2109 // LogStream << "), and " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments() << " line segments, co-incident profiles and their line segments are ";
2110 // for (int nLineSeg = 0; nLineSeg < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments(); nLineSeg++)
2111 // {
2112 // vector<pair<int, int> > prVTmp = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pprVGetPairedCoincidentProfilesForLineSegment(nLineSeg);
2113 // LogStream << "{ ";
2114 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumCoincidentProfilesInLineSegment(nLineSeg); nn++)
2115 // LogStream << prVTmp[nn].first << "[" << prVTmp[nn].second << "] ";
2116 // LogStream << "} ";
2117 // }
2118 // LogStream << endl;
2119 //
2120 // for (int nPoint = 0; nPoint < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize()-1; nPoint++)
2121 // {
2122 // CGeom2DPoint
2123 // Pt1 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint),
2124 // Pt2 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint+1);
2125 //
2126 // if (Pt1 == Pt2)
2127 // LogStream << m_ulIter << ": IDENTICAL POINTS before changes, in profile {" << nThisProfile << "} points = " << nPoint << " and " << nPoint+1 << endl;
2128 // }
2129 // }
2130 // // DEBUG CODE ******************************************************************
2131
2132 return RTN_OK;
2133}
2134
2135//===============================================================================================================================
2137//===============================================================================================================================
2138void CSimulation::TruncateProfileAndAppendNew(int const nCoast, CGeomProfile* pProfileToRetain, int const nMainProfileIntersectLineSeg, vector<CGeom2DPoint> const* pPtVProfileLastPart, vector<vector<pair<int, int>>> const* pprVLineSegLastPart)
2139{
2140 // // DEBUG CODE ****************************************************************
2141 // Get the index numbers of all coincident profiles for the 'main' profile for the line segment in which intersection occurred
2142 // vector<pair<int, int> > prVCoincidentProfilesCHECK1 = *m_VCoast[nCoast].pGetProfile(nMainProfile)->pprVGetPairedCoincidentProfilesForLineSegment(nMainProfileIntersectLineSeg);
2143 // int nNumCoincidentCHECK1 = prVCoincidentProfilesCHECK1.size();
2144 //
2145 // LogStream << "\tTruncating profile {" << nMainProfile << "}, intersection is at [" << dIntersectX << ", " << dIntersectY << "] in line segment " << nMainProfileIntersectLineSeg << endl;
2146 // for (int nn = 0; nn < nNumCoincidentCHECK1; nn++)
2147 // {
2148 // int nThisProfile = prVCoincidentProfilesCHECK1[nn].first;
2149 // LogStream << "\tBEFORE TruncateProfileAndAppendNew(): " << (nThisProfile == nMainProfile ? "MAIN" : "COINCIDENT") << " to-truncate profile {" << nThisProfile << "} has " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize() << " points (";
2150 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize(); nn++)
2151 // LogStream << "[" << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetX() << ", " << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetY() << "] ";
2152 // LogStream << "), and " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments() << " line segments, co-incident profiles are ";
2153 // for (int nLineSeg = 0; nLineSeg < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments(); nLineSeg++)
2154 // {
2155 // vector<pair<int, int> > prVTmp = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pprVGetPairedCoincidentProfilesForLineSegment(nLineSeg);
2156 // LogStream << "{ ";
2157 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumCoincidentProfilesInLineSegment(nLineSeg); nn++)
2158 // LogStream << prVTmp[nn].first << "[" << prVTmp[nn].second << "] ";
2159 // LogStream << "} ";
2160 // }
2161 // LogStream << endl;
2162 //
2163 // for (int nPoint = 0; nPoint < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize()-1; nPoint++)
2164 // {
2165 // CGeom2DPoint
2166 // Pt1 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint),
2167 // Pt2 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint+1);
2168 //
2169 // if (Pt1 == Pt2)
2170 // LogStream << m_ulIter << ": IDENTICAL POINTS before changes, in profile {" << nThisProfile << "} points = " << nPoint << " and " << nPoint+1 << endl;
2171 // }
2172 // }
2173 // LogStream << "\tPart-profile to append is ";
2174 // for (int mm = 0; mm < pPtVProfileLastPart->size(); mm++)
2175 // LogStream << "[" << pPtVProfileLastPart->at(mm).dGetX() << ", " << pPtVProfileLastPart->at(mm).dGetY() << "] ";
2176 // LogStream << endl;
2177 // LogStream << "\tPart line-segment to append is ";
2178 // for (int mm = 0; mm < pprVLineSegLastPart->size(); mm++)
2179 // {
2180 // vector<pair<int, int> > prVTmp = pprVLineSegLastPart->at(mm);
2181 // LogStream << "{ ";
2182 // for (int nn = 0; nn < prVTmp.size(); nn++)
2183 // LogStream << prVTmp[nn].first << "[" << prVTmp[nn].second << "] ";
2184 // LogStream << "} ";
2185 // }
2186 // LogStream << endl;
2187 // // DEBUG CODE ******************************************************************
2188
2189 // Get the index numbers of all coincident profiles for the 'main' profile for the line segment in which intersection occurs
2190 vector<pair<int, int>> prVCoincidentProfiles = *pProfileToRetain->pprVGetPairedCoincidentProfilesForLineSegment(nMainProfileIntersectLineSeg);
2191 int const nNumCoincident = static_cast<int>(prVCoincidentProfiles.size());
2192
2193 for (int nn = 0; nn < nNumCoincident; nn++)
2194 {
2195 // Do this for the main to-truncate profile, and do the same for all its co-incident profiles
2196 int const nThisProfile = prVCoincidentProfiles[nn].first;
2197 int const nThisProfileLineSeg = prVCoincidentProfiles[nn].second;
2198 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfile(nThisProfile);
2199
2200 // if (nThisProfile == nMainProfile)
2201 // assert(nThisProfileLineSeg == nMainProfileIntersectLineSeg);
2202
2203 // Truncate the profile
2204 // LogStream << "\tTruncating " << (nThisProfile == nMainProfile ? "MAIN" : "COINCIDENT") << " to-truncate profile {" << nThisProfile << "} at line segment " << nThisProfileLineSeg+1 << endl;
2205 pThisProfile->TruncateProfile(nThisProfileLineSeg + 1);
2206
2207 // Reduce the number of line segments for this profile
2208 pThisProfile->TruncateLineSegments(nThisProfileLineSeg + 1);
2209
2210 // Append the profile points from the last part of the retain-profile
2211 for (unsigned int mm = 0; mm < pPtVProfileLastPart->size(); mm++)
2212 {
2213 CGeom2DPoint const Pt = pPtVProfileLastPart->at(mm);
2214 pThisProfile->AppendPointInProfile(&Pt);
2215 }
2216
2217 // Append the line segments, and their co-incident profile numbers, from the last part of the retain-profile
2218 for (unsigned int mm = 0; mm < pprVLineSegLastPart->size(); mm++)
2219 {
2220 vector<pair<int, int>> prVTmp = pprVLineSegLastPart->at(mm);
2221
2222 pThisProfile->AppendLineSegment(&prVTmp);
2223 }
2224
2225 // Fix the line seg numbers for this profile
2226 vector<int> nVProf;
2227 vector<int> nVProfsLineSeg;
2228
2229 for (int nSeg = 0; nSeg < pThisProfile->nGetNumLineSegments(); nSeg++)
2230 {
2231 for (int nCoinc = 0; nCoinc < pThisProfile->nGetNumCoincidentProfilesInLineSegment(nSeg); nCoinc++)
2232 {
2233 int const nProf = pThisProfile->nGetProf(nSeg, nCoinc);
2234 int const nProfsLineSeg = pThisProfile->nGetProfsLineSeg(nSeg, nCoinc);
2235
2236 auto it = find(nVProf.begin(), nVProf.end(), nProf);
2237
2238 if (it == nVProf.end())
2239 {
2240 // Not found
2241 nVProf.push_back(nProf);
2242 nVProfsLineSeg.push_back(nProfsLineSeg);
2243 }
2244
2245 else
2246 {
2247 // Found
2248 int const nPos = static_cast<int>(it - nVProf.begin());
2249 int nNewProfsLineSeg = nVProfsLineSeg[nPos];
2250 nNewProfsLineSeg++;
2251
2252 nVProfsLineSeg[nPos] = nNewProfsLineSeg;
2253 pThisProfile->SetProfsLineSeg(nSeg, nCoinc, nNewProfsLineSeg);
2254 }
2255 }
2256 }
2257
2258 // assert(pThisProfile->nGetProfileSize() > 1);
2259 }
2260
2261 // // DEBUG CODE ****************************************************************
2262 // Get the index numbers of all coincident profiles for the 'main' to-truncate profile for the line segment in which intersection occurred
2263 // vector<pair<int, int> > prVToTruncateCoincidentProfilesCHECK2 = *m_VCoast[nCoast].pGetProfile(nMainProfile)->pprVGetPairedCoincidentProfilesForLineSegment(nMainProfileIntersectLineSeg);
2264 // int nNumToTruncateCoincidentCHECK2 = prVToTruncateCoincidentProfilesCHECK2.size();
2265 // for (int nn = 0; nn < nNumToTruncateCoincidentCHECK2; nn++)
2266 // {
2267 // int nThisProfile = prVToTruncateCoincidentProfilesCHECK2[nn].first;
2268 // LogStream << "\tAFTER TruncateProfileAndAppendNew(): " << (nThisProfile == nMainProfile ? "MAIN" : "COINCIDENT") << " to-truncate profile {" << nThisProfile << "} has " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize() << " points (";
2269 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize(); nn++)
2270 // LogStream << "[" << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetX() << ", " << m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nn)->dGetY() << "] ";
2271 // LogStream << "), and " << m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments() << " line segments, co-incident profiles are ";
2272 // for (int mm = 0; mm < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumLineSegments(); mm++)
2273 // {
2274 // vector<pair<int, int> > prVTmp = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pprVGetPairedCoincidentProfilesForLineSegment(mm);
2275 // LogStream << "{ ";
2276 // for (int nn = 0; nn < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetNumCoincidentProfilesInLineSegment(mm); nn++)
2277 // LogStream << prVTmp[nn].first << "[" << prVTmp[nn].second << "] ";
2278 // LogStream << "} ";
2279 // }
2280 // LogStream << endl;
2281 //
2282 // for (int nPoint = 0; nPoint < m_VCoast[nCoast].pGetProfile(nThisProfile)->nGetProfileSize()-1; nPoint++)
2283 // {
2284 // CGeom2DPoint
2285 // Pt1 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint),
2286 // Pt2 = *m_VCoast[nCoast].pGetProfile(nThisProfile)->pPtGetPointInProfile(nPoint+1);
2287 //
2288 // if (Pt1 == Pt2)
2289 // LogStream << m_ulIter << ": IDENTICAL POINTS before changes, in profile {" << nThisProfile << "} points = " << nPoint << " and " << nPoint+1 << endl;
2290 // }
2291 // }
2292 // // DEBUG CODE ******************************************************************
2293}
2294
Contains CGeom2DPoint definitions.
Contains CGeom2DIPoint definitions.
vector< CGeom2DPoint > * pPtVGetPoints(void)
Returns the address of the vector which represents this 2D shape.
Definition 2d_shape.cpp:130
Geometry class used to represent 2D point objects with integer coordinates.
Definition 2di_point.h:29
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
Definition 2di_point.cpp:55
void SetXY(int const, int const)
The two integer parameters set values for the CGeom2DIPoint object's X and Y coordinates.
Definition 2di_point.cpp:85
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Definition 2di_point.cpp:49
Geometry class used to represent 2D point objects with floating-point coordinates.
Definition 2d_point.h:27
void SetY(double const)
The double parameter sets a value for the CGeom2DIPoint object's Y coordinate.
Definition 2d_point.cpp:65
double dGetY(void) const
Returns the CGeom2DPoint object's double Y coordinate.
Definition 2d_point.cpp:53
double dGetX(void) const
Returns the CGeom2DPoint object's double X coordinate.
Definition 2d_point.cpp:47
void SetX(double const)
The double parameter sets a value for the CGeom2DIPoint object's X coordinate.
Definition 2d_point.cpp:59
int nGetProfsLineSeg(int const, int const) const
Returns the profile's own line segment, given a line segment and the index of the co-incident profile...
bool bFindProfileInCoincidentProfilesOfLastLineSegment(int const)
Returns true if the given profile number is amongst the coincident profiles of the CGeomMultiLine obj...
int nGetNumCoincidentProfilesInLineSegment(int const)
Returns the count of coincident profiles in a specified line segment, or -1 if the line segment does ...
void AppendLineSegment(void)
Appends a new empty line segment.
int nGetProf(int const, int const) const
Returns the profile number, given a line segment and the index of the co-incident profile for that li...
void TruncateLineSegments(int const)
Cuts short the number of line segments.
vector< vector< pair< int, int > > > prVVGetAllLineSegAfter(int const)
Returns a vector of the line segments which succeed the specified line segment number.
void AppendCoincidentProfileToLineSegments(pair< int, int > const)
Appends a coincident profile pair to the CGeomMultiLine object's final line segment.
void SetProfsLineSeg(int const, int const, int const)
Sets a profile's own line segment number, given a line segment and the index of the co-incident profi...
vector< pair< int, int > > * pprVGetPairedCoincidentProfilesForLineSegment(int const)
Returns a vector of pairs (a line segment)
void AddCoincidentProfileToExistingLineSegment(int const, int const, int const)
Adds a coincident profile to a pre-existing line segment of the CGeomMultiLine object.
int nGetNumLineSegments(void) const
Appends a line segment which then inherits from the preceding line segments.
Geometry class used to represent coast profile objects.
Definition profile.h:37
void TruncateProfile(int const)
Truncates the profile's CGeomLine (external CRS points)
Definition profile.cpp:341
void SetEndOfCoast(bool const)
Sets a switch to indicate whether this is an end-of-coast profile.
Definition profile.cpp:122
void AppendPointInProfile(double const, double const)
Appends a point (external CRS) to the profile.
Definition profile.cpp:310
int nGetProfileID(void) const
Returns the profile's this-coast ID.
Definition profile.cpp:80
void SetProfileDeepWaterWavePeriod(double const)
Sets the deep-water wave period for this profile.
Definition profile.cpp:591
bool bIsPointInProfile(double const, double const)
Removes a line segment from the profile.
Definition profile.cpp:386
void SetHitCoast(bool const)
Sets a switch which indicates whether this profile has hit a coast.
Definition profile.cpp:170
bool bProfileOKIncTruncated(void) const
Returns true if this is a problem-free profile, and is not a start-of-coast and is not an end-of-coas...
Definition profile.cpp:248
void SetStartOfCoast(bool const)
Sets a switch to indicate whether this is a start-of-coast profile.
Definition profile.cpp:110
void SetUpCoastAdjacentProfile(CGeomProfile *)
Sets the up-coast adjacent profile.
Definition profile.cpp:473
void SetHitAnotherProfile(bool const)
Sets a switch which indicates whether this profile hits another profile badly.
Definition profile.cpp:218
vector< CGeom2DPoint > PtVGetThisPointAndAllAfter(int const)
Returns a given external CRS point from the profile, and all points after this.
Definition profile.cpp:373
bool bInsertIntersection(double const, double const, int const)
Inserts an intersection (at a point specified in external CRS, with a line segment) into the profile.
Definition profile.cpp:322
void SetTruncatedSameCoast(bool const)
Sets a switch which indicates whether this profile is truncated, due to hitting another profile from ...
Definition profile.cpp:194
bool bProfileOK(void) const
Returns true if this is a problem-free profile, and is not a start-of-coast and is not an end-of-coas...
Definition profile.cpp:230
void SetProfileDeepWaterWaveHeight(double const)
Sets the deep-water wave height for this profile.
Definition profile.cpp:567
void SetDownCoastAdjacentProfile(CGeomProfile *)
Sets the down-coast adjacent profile.
Definition profile.cpp:485
bool bIsIntervention(void) const
Returns true if this is an intervention profile.
Definition profile.cpp:603
void SetHitIntervention(bool const)
Sets a switch which indicates whether this profile has hit an intervention.
Definition profile.cpp:158
CGeom2DPoint * pPtGetPointInProfile(int const)
Returns a single point (external CRS) from the profile.
Definition profile.cpp:367
void AppendCellInProfile(CGeom2DIPoint const *)
Appends a cell (grid CRS) to the profile.
Definition profile.cpp:497
void SetPointInProfile(int const, double const, double const)
Sets a single point (external CRS) in the profile.
Definition profile.cpp:303
int nGetProfileSize(void) const
Returns the number of external CRS points in the profile (only two, initally; and always just two for...
Definition profile.cpp:361
CGeom2DIPoint * pPtiGetStartPoint(void)
Returns a pointer to the location of the cell (grid CRS) on which the profile starts.
Definition profile.cpp:92
bool bStartOfCoast(void) const
Returns the switch to indicate whether this is a start-of-coast profile.
Definition profile.cpp:116
void SetEndPoint(CGeom2DIPoint const *)
Sets the the location of the cell (grid CRS) on which the profile ends.
Definition profile.cpp:98
void SetHitLand(bool const)
Sets a switch which indicates whether this profile has hit land.
Definition profile.cpp:146
void SetPointsInProfile(vector< CGeom2DPoint > const *)
Sets points (external CRS) in the profile. Note that only two points, the start and end point,...
Definition profile.cpp:297
bool bEndOfCoast(void) const
Returns the switch to indicate whether this is an end-of-coast profile.
Definition profile.cpp:128
void SetProfileDeepWaterWaveAngle(double const)
Sets the deep-water wave orientation for this profile.
Definition profile.cpp:579
void SetTooShort(bool const)
Sets a switch which indicates whether this profile is too short to be useful.
Definition profile.cpp:182
int m_nLogFileDetail
The level of detail in the log file output. Can be LOG_FILE_LOW_DETAIL, LOG_FILE_MIDDLE_DETAIL,...
Definition simulation.h:574
double m_dThisIterSWL
The still water level for this timestep (this includes tidal changes and any long-term SWL change)
Definition simulation.h:712
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
int m_nXGridSize
The size of the grid in the x direction.
Definition simulation.h:457
void LocateAndCreateProfiles(int const, int &, vector< bool > *, vector< pair< int, double > > const *)
For a single coastline, locate the start points for all coastline-normal profiles (except the grid-ed...
CGeom2DIPoint PtiExtCRSToGridRound(CGeom2DPoint const *) const
Transforms a pointer to a CGeom2DPoint in the external CRS to the equivalent CGeom2DIPoint in the ras...
ofstream LogStream
vector< CRWCoast > m_VCoast
The coastline objects.
double m_dCoastNormalLength
Length of the coastline-normal profiles, in m.
Definition simulation.h:826
void MarkProfilesOnGrid(int const, int &)
For this coastline, marks all coastline-normal profiles (apart from the two 'special' ones at the sta...
static CGeom2DPoint PtChooseEndPoint(int const, CGeom2DPoint const *, CGeom2DPoint const *, double const, double const, double const, double const)
Choose which end point to use for the coastline-normal profile.
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:460
normal_distribution< double > m_dGetFromUnitNormalDist
c++11 unit normal distribution (mean = 0, stdev = 1)
void CreateRasterizedProfile(int const, CGeomProfile *, vector< CGeom2DIPoint > *, vector< bool > *, bool &, bool &, bool &, bool &, bool &, bool &)
Given a pointer to a coastline-normal profile, returns an output vector of cells which are 'under' ev...
static CGeom2DPoint PtAverage(CGeom2DPoint const *, CGeom2DPoint const *)
Returns a point (external CRS) which is the average of (i.e. is midway between) two other external CR...
void KeepWithinValidGrid(int &, int &) const
int m_nCoastNormalSpacing
Average spacing between coastline normals, measured in cells.
Definition simulation.h:490
void MergeProfilesAtFinalLineSegments(int const, CGeomProfile *, CGeomProfile *, int const, int const, double const, double const, double const, double const)
Merges two profiles which intersect at their final (most seaward) line segments, seaward of their poi...
int m_nCoastNormalInterventionSpacing
Average spacing between coastline normals on interventions, measured in cells.
Definition simulation.h:493
void CheckForIntersectingProfiles(void)
Checks all coastline-normal profiles for intersection, and modifies those that intersect.
double dGridCentroidYToExtCRSY(int const) const
Given the integer Y-axis ordinate of a cell in the raster grid CRS, returns the external CRS Y-axis o...
Definition gis_utils.cpp:78
int nCreateProfile(int const, int const, int const, int const, bool const, CGeom2DIPoint const *)
Creates a single coastline-normal profile (which may be an intervention profile)
vector< CGeom2DIPoint > m_VEdgeCell
Edge cells.
static bool bCheckForIntersection(CGeomProfile *const, CGeomProfile *const, int &, int &, double &, double &, double &, double &)
Checks all line segments of a pair of coastline-normal profiles for intersection. If the lines inters...
int nGetCoastNormalEndPoint(int const, int const, int const, CGeom2DPoint const *, double const, CGeom2DPoint *, CGeom2DIPoint *, bool const)
Finds the end point of a coastline-normal line, given the start point on the vector coastline....
void TruncateProfileAndAppendNew(int const, CGeomProfile *, int const, vector< CGeom2DPoint > const *, vector< vector< pair< int, int > > > const *)
Truncate a profile at the point of intersection, and do the same for all its co-incident profiles.
double dExtCRSXToGridX(double const) const
int nInsertPointIntoProfilesIfNeededThenUpdate(int const, CGeomProfile *, double const, double const, int const, CGeomProfile *, int const, bool const)
Inserts an intersection point into the profile that is to be retained, if that point is not already p...
default_random_engine m_Rand[NUMBER_OF_RNGS]
The c++11 random number generators.
int nCreateAllProfiles(void)
Create coastline-normal profiles for all coastlines. The first profiles are created 'around' the most...
double dExtCRSYToGridY(double const) const
Transforms a Y-axis ordinate in the external CRS to the equivalent Y-axis ordinate in the raster grid...
double m_dCoastNormalRandSpacingFactor
Random factor for spacing of along-coast normals.
Definition simulation.h:952
bool bIsWithinValidGrid(int const, int const) const
static void AppendEnsureNoGap(vector< CGeom2DIPoint > *, CGeom2DIPoint const *)
Appends a CGeom2DIPoint to a vector<CGeom2DIPoint>, making sure that the new end point touches the pr...
Definition utils.cpp:2526
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:598
double dGridCentroidXToExtCRSX(int const) const
Definition gis_utils.cpp:68
int nLocateAndCreateGridEdgeProfile(bool const, int const, int &)
Creates a 'special' profile at each end of a coastline, at the edge of the raster grid....
void TruncateOneProfileRetainOtherProfile(int const, CGeomProfile *, CGeomProfile *, double, double, int, int, bool const)
Truncates one intersecting profile at the point of intersection, and retains the other profile.
double m_dCellSide
Length of a cell side (in external CRS units)
Definition simulation.h:658
double m_dDepthOfClosure
Depth of closure (in m) TODO 007 can be calculated using Hallermeier, R.J. (1978) or Birkemeier (1985...
Definition simulation.h:817
int nCheckAndMarkAllProfiles(void)
Check all coastline-normal profiles and modify the profiles if they intersect, then mark valid profil...
vector< int > m_VEdgeCellEdge
The grid edge that each edge cell belongs to.
This file contains global definitions for CoastalME.
int const INT_NODATA
Definition cme.h:476
T tMin(T a, T b)
Definition cme.h:1265
double const TOLERANCE
Definition cme.h:823
int const RTN_ERR_COAST_CANT_FIND_EDGE_CELL
Definition cme.h:750
string const ERR
Definition cme.h:902
int const DIRECTION_DOWNCOAST
Definition cme.h:506
bool const ACCEPT_TRUNCATED_PROFILES
Definition cme.h:457
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:491
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1303
T tMax(T a, T b)
Definition cme.h:1252
int const RTN_ERR_PROFILE_END_INSUFFICIENT_DEPTH
Definition cme.h:722
int const PROFILE_CHECK_DIST_FROM_COAST
Definition cme.h:799
string const WARN
Definition cme.h:903
int const RTN_ERR_NO_PROFILES_2
Definition cme.h:724
int const RTN_ERR_PROFILE_ENDPOINT_IS_INLAND
Definition cme.h:720
int const LOG_FILE_ALL
Definition cme.h:493
int const RTN_OK
Definition cme.h:694
int const RTN_ERR_NO_PROFILES_1
Definition cme.h:723
int const RTN_ERR_CANNOT_INSERT_POINT
Definition cme.h:741
int const DIRECTION_UPCOAST
Definition cme.h:507
string strDblToStr(const T &t)
Definition cme.h:1291
int const RTN_ERR_NO_SOLUTION_FOR_ENDPOINT
Definition cme.h:721
int const RIGHT_HANDED
Definition cme.h:511
T tAbs(T a)
Definition cme.h:1277
int const LF_CAT_INTERVENTION
Definition cme.h:545
int const LEFT_HANDED
Definition cme.h:512
Contains CRWCoast definitions.
Contains CSimulation definitions.
int nRound(double const d)
Version of the above that returns an int.