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