CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
locate_coast.cpp
Go to the documentation of this file.
1
12
13/*==============================================================================================================================
14
15This file is part of CoastalME, the Coastal Modelling Environment.
16
17CoastalME 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
19This 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
21You 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#include <cfloat>
26
27#include <iostream>
28using std::cerr;
29using std::cout;
30using std::endl;
31using std::ios;
32
33#include <iomanip>
34using std::resetiosflags;
35using std::setiosflags;
36using std::setprecision;
37using std::setw;
38
39#include <string>
40using std::to_string;
41
42#include <stack>
43using std::stack;
44
45#include "cme.h"
46#include "i_line.h"
47#include "line.h"
48#include "simulation.h"
49#include "raster_grid.h"
50#include "coast.h"
51
52//===============================================================================================================================
54//===============================================================================================================================
55int CSimulation::nLocateSeaAndCoasts(int &nValidCoast)
56{
57 // Find all connected sea cells
58 FindAllSeaCells();
59
60 // Find every coastline on the raster grid, mark raster cells, then create the vector coastline
61 int nRet = nTraceAllCoasts(nValidCoast);
62 if (nRet != RTN_OK)
63 return nRet;
64
65 // Have we created any coasts?
66 if (m_VCoast.empty())
67 {
68 cerr << m_ulIter << ": " << ERR << "no coastline located: this iteration SWL = " << m_dThisIterSWL << ", maximum DEM top surface elevation = " << m_dThisIterTopElevMax << ", minimum DEM top surface elevation = " << m_dThisIterTopElevMin << endl;
69 return RTN_ERR_NOCOAST;
70 }
71
72 return RTN_OK;
73}
74
75//===============================================================================================================================
77//===============================================================================================================================
78void CSimulation::FindAllSeaCells(void)
79{
80 // Go along the list of edge cells
81 for (unsigned int n = 0; n < m_VEdgeCell.size(); n++)
82 {
83 if (m_bOmitSearchNorthEdge && m_VEdgeCellEdge[n] == NORTH)
84 continue;
85
86 if (m_bOmitSearchSouthEdge && m_VEdgeCellEdge[n] == SOUTH)
87 continue;
88
89 if (m_bOmitSearchWestEdge && m_VEdgeCellEdge[n] == WEST)
90 continue;
91
92 if (m_bOmitSearchEastEdge && m_VEdgeCellEdge[n] == EAST)
93 continue;
94
95 int nX = m_VEdgeCell[n].nGetX();
96 int nY = m_VEdgeCell[n].nGetY();
97
98 // if ((m_pRasterGrid->m_Cell[nX][nY].bIsInundated()) && (m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth() == 0))
99 if ((m_pRasterGrid->m_Cell[nX][nY].bIsInundated()) && (bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth(), 0.0, TOLERANCE)))
100
101 // This edge cell is below SWL and sea depth remains set to zero
102 FloodFillSea(nX, nY);
103 }
104}
105
106//===============================================================================================================================
108//===============================================================================================================================
109void CSimulation::FloodFillSea(int const nXStart, int const nYStart)
110{
111 // For safety check
112 int nRoundLoopMax = m_nXGridSize * m_nYGridSize;
113
114 // Create an empty stack
115 stack<CGeom2DIPoint> PtiStack;
116
117 // Start at the given edge cell, push this onto the stack
118 PtiStack.push(CGeom2DIPoint(nXStart, nYStart));
119
120 // Then do the flood fill loop until there are no more cell coordinates on the stack
121 int nRoundLoop = 0;
122 while (! PtiStack.empty())
123 {
124 // Safety check
125 if (nRoundLoop++ > nRoundLoopMax)
126 break;
127
128 CGeom2DIPoint Pti = PtiStack.top();
129 PtiStack.pop();
130
131 int nX = Pti.nGetX();
132 int nY = Pti.nGetY();
133
134 while ((nX >= 0) && (! m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue()) && (m_pRasterGrid->m_Cell[nX][nY].bIsInundated()))
135 nX--;
136
137 nX++;
138
139 bool bSpanAbove = false;
140 bool bSpanBelow = false;
141
142 while ((nX < m_nXGridSize) && (! m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue()) && (m_pRasterGrid->m_Cell[nX][nY].bIsInundated()) && (bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth(), 0.0, TOLERANCE)))
143 {
144 // Set the sea depth for this cell
145 m_pRasterGrid->m_Cell[nX][nY].SetSeaDepth();
146
147 CRWCellLandform* pLandform = m_pRasterGrid->m_Cell[nX][nY].pGetLandform();
148 int nCat = pLandform->nGetLFCategory();
149
150 // Have we had sediment input here?
151 if ((nCat == LF_CAT_SEDIMENT_INPUT) || (nCat == LF_CAT_SEDIMENT_INPUT_SUBMERGED) || (nCat == LF_CAT_SEDIMENT_INPUT_NOT_SUBMERGED))
152 {
153 if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
154 {
155 pLandform->SetLFCategory(LF_CAT_SEDIMENT_INPUT_SUBMERGED);
156 m_pRasterGrid->m_Cell[nX][nY].SetInContiguousSea();
157
158 // Set this sea cell to have deep water (off-shore) wave orientation and height, will change this later for cells closer to the shoreline if we have on-shore waves
159 m_pRasterGrid->m_Cell[nX][nY].SetWaveValuesToDeepWaterWaveValues();
160 }
161 else
162 {
163 pLandform->SetLFCategory(LF_CAT_SEDIMENT_INPUT_NOT_SUBMERGED);
164 }
165 }
166 else
167 {
168 // No sediment input here, just mark as sea
169 m_pRasterGrid->m_Cell[nX][nY].SetInContiguousSea();
170 pLandform->SetLFCategory(LF_CAT_SEA);
171
172 // Set this sea cell to have deep water (off-shore) wave orientation and height, will change this later for cells closer to the shoreline if we have on-shore waves
173 m_pRasterGrid->m_Cell[nX][nY].SetWaveValuesToDeepWaterWaveValues();
174 }
175
176 // Now sort out the x-y extremities of the contiguous sea for the bounding box (used later in wave propagation)
177 if (nX < m_nXMinBoundingBox)
178 m_nXMinBoundingBox = nX;
179
180 if (nX > m_nXMaxBoundingBox)
181 m_nXMaxBoundingBox = nX;
182
183 if (nY < m_nYMinBoundingBox)
184 m_nYMinBoundingBox = nY;
185
186 if (nY > m_nYMaxBoundingBox)
187 m_nYMaxBoundingBox = nY;
188
189 // Update count
190 m_ulThisIterNumSeaCells++;
191
192 if ((! bSpanAbove) && (nY > 0) && (! m_pRasterGrid->m_Cell[nX][nY-1].bBasementElevIsMissingValue()) && (m_pRasterGrid->m_Cell[nX][nY-1].bIsInundated()))
193 {
194 PtiStack.push(CGeom2DIPoint(nX, nY-1));
195 bSpanAbove = true;
196 }
197 else if (bSpanAbove && (nY > 0) && (! m_pRasterGrid->m_Cell[nX][nY-1].bBasementElevIsMissingValue()) && (! m_pRasterGrid->m_Cell[nX][nY-1].bIsInundated()))
198 {
199 bSpanAbove = false;
200 }
201
202 if ((! bSpanBelow) && (nY < m_nYGridSize-1) && (! m_pRasterGrid->m_Cell[nX][nY+1].bBasementElevIsMissingValue()) && (m_pRasterGrid->m_Cell[nX][nY+1].bIsInundated()))
203 {
204 PtiStack.push(CGeom2DIPoint(nX, nY+1));
205 bSpanBelow = true;
206 }
207 else if (bSpanBelow && (nY < m_nYGridSize-1) && (! m_pRasterGrid->m_Cell[nX][nY+1].bBasementElevIsMissingValue()) && (! m_pRasterGrid->m_Cell[nX][nY+1].bIsInundated()))
208 {
209 bSpanBelow = false;
210 }
211
212 nX++;
213 }
214 }
215
216 // // DEBUG CODE ===========================================================================================================
217 // string strOutFile = m_strOutPath + "is_contiguos_sea_";
218 // strOutFile += to_string(m_ulIter);
219 // strOutFile += ".tif";
220
221 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
222 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
223 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
224 // pDataSet->SetGeoTransform(m_dGeoTransform);
225 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
226 // int n = 0;
227 // for (int nY = 0; nY < m_nYGridSize; nY++)
228 // {
229 // for (int nX = 0; nX < m_nXGridSize; nX++)
230 // {
231 // pdRaster[n++] = m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea();
232 // }
233 // }
234
235 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
236 // pBand->SetNoDataValue(m_dMissingValue);
237 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
238 // if (nRet == CE_Failure)
239 // return;
240
241 // GDALClose(pDataSet);
242 // delete[] pdRaster;
243
244 // string strOutFile = m_strOutPath + "is_inundated_";
245 // strOutFile += to_string(m_ulIter);
246 // strOutFile += ".tif";
247
248 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
249 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
250 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
251 // pDataSet->SetGeoTransform(m_dGeoTransform);
252 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
253
254 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
255 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
256 // pDataSet->SetGeoTransform(m_dGeoTransform);
257
258 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
259 // int n = 0;
260 // for (int nY = 0; nY < m_nYGridSize; nY++)
261 // {
262 // for (int nX = 0; nX < m_nXGridSize; nX++)
263 // {
264 // pdRaster[n++] = m_pRasterGrid->m_Cell[nX][nY].bIsInundated();
265 // }
266 // }
267
268 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
269 // pBand = pDataSet->GetRasterBand(1);
270 // pBand->SetNoDataValue(m_dMissingValue);
271 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
272 // if (nRet == CE_Failure)
273 // return;
274
275 // GDALClose(pDataSet);
276 // delete[] pdRaster;
277 // // DEBUG CODE ===========================================================================================================
278
279 // LogStream << m_ulIter << ": flood fill of sea from [" << nXStart << "][" << nYStart << "] = {" << dGridCentroidXToExtCRSX(nXStart) << ", " << dGridCentroidYToExtCRSY(nYStart) << "} with SWL = " << m_dThisIterSWL << ", " << m_ulThisIterNumSeaCells << " of " << m_ulNumCells << " cells now marked as sea (" << fixed << setprecision(3) << 100.0 * m_ulThisIterNumSeaCells / m_ulNumCells << " %)" << endl;
280
281 // LogStream << " m_nXMinBoundingBox = " << m_nXMinBoundingBox << " m_nXMaxBoundingBox = " << m_nXMaxBoundingBox << " m_nYMinBoundingBox = " << m_nYMinBoundingBox << " m_nYMaxBoundingBox = " << m_nYMaxBoundingBox << endl;
282}
283
284//===============================================================================================================================
286//===============================================================================================================================
287int CSimulation::nTraceAllCoasts(int &nValidCoast)
288{
289 if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
290 LogStream << m_ulIter << ": Tracing coast" << endl;
291
292 vector<bool> VbPossibleStartCellLHEdge;
293 vector<bool> VbTraced;
294 vector<int> VnSearchDirection;
295 vector<CGeom2DIPoint> V2DIPossibleStartCell;
296
297 // Go along the list of edge cells and look for possible coastline start cells
298 for (unsigned int n = 0; n < m_VEdgeCell.size() - 1; n++)
299 {
300 if (m_bOmitSearchNorthEdge && (m_VEdgeCellEdge[n] == NORTH || m_VEdgeCellEdge[n + 1] == NORTH))
301 continue;
302
303 if (m_bOmitSearchSouthEdge && (m_VEdgeCellEdge[n] == SOUTH || m_VEdgeCellEdge[n + 1] == SOUTH))
304 continue;
305
306 if (m_bOmitSearchWestEdge && (m_VEdgeCellEdge[n] == WEST || m_VEdgeCellEdge[n + 1] == WEST))
307 continue;
308
309 if (m_bOmitSearchEastEdge && (m_VEdgeCellEdge[n] == EAST || m_VEdgeCellEdge[n + 1] == EAST))
310 continue;
311
312 int nXThis = m_VEdgeCell[n].nGetX();
313 int nYThis = m_VEdgeCell[n].nGetY();
314 int nXNext = m_VEdgeCell[n + 1].nGetX();
315 int nYNext = m_VEdgeCell[n + 1].nGetY();
316
317 // Get "Is it sea?" information for 'this' and 'next' cells
318 bool bThisCellIsSea = m_pRasterGrid->m_Cell[nXThis][nYThis].bIsInContiguousSea();
319 bool bNextCellIsSea = m_pRasterGrid->m_Cell[nXNext][nYNext].bIsInContiguousSea();
320
321 // Are we at a coast?
322 if ((! bThisCellIsSea) && bNextCellIsSea)
323 {
324 // 'This' cell is just inland, has it already been flagged as a possible start for a coastline (even if this subsequently 'failed' as a coastline)?
325 if (! m_pRasterGrid->m_Cell[nXThis][nYThis].bIsPossibleCoastStartCell())
326 {
327 // It has not, so flag it
328 m_pRasterGrid->m_Cell[nXThis][nYThis].SetPossibleCoastStartCell();
329 if (m_nLogFileDetail >= LOG_FILE_ALL)
330 LogStream << m_ulIter << ": Flagging [" << nXThis << "][" << nYThis << "] = {" << dGridCentroidXToExtCRSX(nXThis) << ", " << dGridCentroidYToExtCRSY(nYThis) << "} as possible coast start cell, with left_handed edge" << endl;
331
332 // And save it
333 V2DIPossibleStartCell.push_back(CGeom2DIPoint(nXThis, nYThis));
334 VbPossibleStartCellLHEdge.push_back(true);
335 VnSearchDirection.push_back(nGetOppositeDirection(m_VEdgeCellEdge[n]));
336 VbTraced.push_back(false);
337 }
338 }
339 else if (bThisCellIsSea && (! bNextCellIsSea))
340 {
341 // The 'next' cell is just inland, has it already been flagged as a possible start for a coastline (even if this subsequently 'failed' as a coastline)?
342 if (! m_pRasterGrid->m_Cell[nXNext][nYNext].bIsPossibleCoastStartCell())
343 {
344 // It has not, so flag it
345 m_pRasterGrid->m_Cell[nXNext][nYNext].SetPossibleCoastStartCell();
346 if (m_nLogFileDetail >= LOG_FILE_ALL)
347 LogStream << m_ulIter << ": Flagging [" << nXNext << "][" << nYNext << "] = {" << dGridCentroidXToExtCRSX(nXNext) << ", " << dGridCentroidYToExtCRSY(nYNext) << "} as possible coast start cell, with right_handed edge" << endl;
348
349 // And save it
350 V2DIPossibleStartCell.push_back(CGeom2DIPoint(nXNext, nYNext));
351 VbPossibleStartCellLHEdge.push_back(false);
352 VnSearchDirection.push_back(nGetOppositeDirection(m_VEdgeCellEdge[n + 1]));
353 VbTraced.push_back(false);
354 }
355 }
356 }
357
358 for (unsigned int n = 0; n < V2DIPossibleStartCell.size(); n++)
359 {
360 if (! VbTraced[n])
361 {
362 int nRet = 0;
363 if (VbPossibleStartCellLHEdge[n])
364 {
365 nRet = nTraceCoastLine(n, VnSearchDirection[n], LEFT_HANDED, &VbTraced, &V2DIPossibleStartCell);
366 }
367 else
368 {
369 nRet = nTraceCoastLine(n, VnSearchDirection[n], RIGHT_HANDED, &VbTraced, &V2DIPossibleStartCell);
370 }
371
372 if (nRet == RTN_OK)
373 {
374 // We have a valid coastline starting from this possible start cell
375 VbTraced[n] = true;
376 nValidCoast++;
377 }
378 }
379 }
380
381 if (nValidCoast > 0)
382 return RTN_OK;
383 else
384 {
385 cerr << m_ulIter << ": no valid coasts found" << endl;
386 return RTN_ERR_TRACING_COAST;
387 }
388}
389
390//===============================================================================================================================
392//===============================================================================================================================
393int CSimulation::nTraceCoastLine(unsigned int const nTraceFromStartCellIndex, int const nStartSearchDirection, int const nHandedness, vector<bool>* pVbTraced, vector<CGeom2DIPoint> const* pV2DIPossibleStartCell)
394{
395 bool bHitStartCell = false;
396 bool bAtCoast = false;
397 bool bHasLeftStartEdge = false;
398 bool bTooLong = false;
399 bool bOffEdge = false;
400 bool bRepeating = false;
401
402 int nStartX = pV2DIPossibleStartCell->at(nTraceFromStartCellIndex).nGetX();
403 int nStartY = pV2DIPossibleStartCell->at(nTraceFromStartCellIndex).nGetY();
404 int nX = nStartX;
405 int nY = nStartY;
406 int nSearchDirection = nStartSearchDirection;
407 int nRoundLoop = -1;
408 // nThisLen = 0;
409 // nLastLen = 0,
410 // nPreLastLen = 0;
411
412 // Temporary coastline as integer points (grid CRS)
413 CGeomILine ILTempGridCRS;
414
415 // Mark the start cell as coast and add it to the vector object
416 m_pRasterGrid->m_Cell[nStartX][nStartY].SetAsCoastline(true);
417 CGeom2DIPoint PtiStart(nStartX, nStartY);
418 ILTempGridCRS.Append(&PtiStart);
419
420 // Start at this grid-edge point and trace the rest of the coastline using the 'wall follower' rule for maze traversal, trying to keep next to cells flagged as sea
421 do
422 {
423 // // DEBUG CODE ==============================================================================================================
424 // LogStream << "Now at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
425 // LogStream << "ILTempGridCRS is now:" << endl;
426 // for (int n = 0; n < ILTempGridCRS.nGetSize(); n++)
427 // LogStream << "[" << ILTempGridCRS[n].nGetX() << "][" << ILTempGridCRS[n].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[n].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[n].nGetY()) << "}" << endl;
428 // LogStream << "=================" << endl;
429 // // DEBUG CODE ==============================================================================================================
430
431 // Safety check
432 if (++nRoundLoop > m_nCoastMax)
433 {
434 bTooLong = true;
435
436 // LogStream << m_ulIter << ": abandoning coastline tracing from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "}, exceeded maximum search length (" << m_nCoastMax << ")" << endl;
437
438 // for (int n = 0; n < ILTempGridCRS.nGetSize(); n++)
439 // LogStream << "[" << ILTempGridCRS[n].nGetX() << "][" << ILTempGridCRS[n].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[n].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[n].nGetY()) << "}" << endl;
440 // LogStream << endl;
441
442 break;
443 }
444
445 // Another safety check
446 if ((nRoundLoop > 10) && (ILTempGridCRS.nGetSize() < 2))
447 {
448 // We've been 10 times round the loop but the coast is still less than 2 coastline points in length, so we must be repeating
449 bRepeating = true;
450
451 break;
452 }
453
454 // OK so far: so have we left the start edge?
455 if (! bHasLeftStartEdge)
456 {
457 // We have not yet left the start edge
458 if (((nStartSearchDirection == SOUTH) && (nY > nStartY)) || ((nStartSearchDirection == NORTH) && (nY < nStartY)) ||
459 ((nStartSearchDirection == EAST) && (nX > nStartX)) || ((nStartSearchDirection == WEST) && (nX < nStartX)))
460 bHasLeftStartEdge = true;
461
462 // Flag this cell to ensure that it is not chosen as a coastline start cell later
463 m_pRasterGrid->m_Cell[nX][nY].SetPossibleCoastStartCell();
464 // LogStream << "Flagging [" << nX << "][" << nY << "] as possible coast start cell NOT YET LEFT EDGE" << endl;
465 }
466
467 // Leave the loop if the vector coastline has left the start edge, then we find a coast cell which is a possible start cell from which a coastline has not yet been traced
468 // if (bHasLeftStartEdge && bAtCoast)
469 {
470 for (unsigned int nn = 0; nn < pVbTraced->size(); nn++)
471 {
472 if ((nn != nTraceFromStartCellIndex) && (!pVbTraced->at(nn)))
473 {
474 // LogStream << "[" << pV2DIPossibleStartCell->at(nn).nGetX() << "][" << pV2DIPossibleStartCell->at(nn).nGetY() << "]" << endl;
475
476 if (bAtCoast && (nX == pV2DIPossibleStartCell->at(nn).nGetX()) && (nY == pV2DIPossibleStartCell->at(nn).nGetY()))
477 {
478 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
479 LogStream << m_ulIter << ": Valid coastline found, traced from [" << nStartX << "][" << nStartY << "] and hit another start cell at [" << nX << "][" << nY << "]" << endl;
480
481 pVbTraced->at(nn) = true;
482 bHitStartCell = true;
483 break;
484 }
485 }
486 }
487 // LogStream << endl;
488 }
489
490 if (bHitStartCell)
491 break;
492
493 // OK now sort out the next iteration of the search
494 int nXSeaward = 0;
495 int nYSeaward = 0;
496 int nSeawardNewDirection = 0;
497 int nXStraightOn = 0;
498 int nYStraightOn = 0;
499 int nXAntiSeaward = 0;
500 int nYAntiSeaward = 0;
501 int nAntiSeawardNewDirection = 0;
502 int nXGoBack = 0;
503 int nYGoBack = 0;
504 int nGoBackNewDirection = 0;
505
506 CGeom2DIPoint Pti(nX, nY);
507
508 // Set up the variables
509 switch (nHandedness)
510 {
511 case RIGHT_HANDED:
512 // The sea is to the right-hand side of the coast as we traverse it. We are just inland, so we need to keep heading right to find the sea
513 switch (nSearchDirection)
514 {
515 case NORTH:
516 // The sea is towards the RHS (E) of the coast, so first try to go right (to the E)
517 nXSeaward = nX + 1;
518 nYSeaward = nY;
519 nSeawardNewDirection = EAST;
520
521 // If can't do this, try to go straight on (to the N)
522 nXStraightOn = nX;
523 nYStraightOn = nY - 1;
524
525 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (W)
526 nXAntiSeaward = nX - 1;
527 nYAntiSeaward = nY;
528 nAntiSeawardNewDirection = WEST;
529
530 // As a last resort, go back (to the S)
531 nXGoBack = nX;
532 nYGoBack = nY + 1;
533 nGoBackNewDirection = SOUTH;
534
535 break;
536
537 case EAST:
538 // The sea is towards the RHS (S) of the coast, so first try to go right (to the S)
539 nXSeaward = nX;
540 nYSeaward = nY + 1;
541 nSeawardNewDirection = SOUTH;
542
543 // If can't do this, try to go straight on (to the E)
544 nXStraightOn = nX + 1;
545 nYStraightOn = nY;
546
547 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (N)
548 nXAntiSeaward = nX;
549 nYAntiSeaward = nY - 1;
550 nAntiSeawardNewDirection = NORTH;
551
552 // As a last resort, go back (to the W)
553 nXGoBack = nX - 1;
554 nYGoBack = nY;
555 nGoBackNewDirection = WEST;
556
557 break;
558
559 case SOUTH:
560 // The sea is towards the RHS (W) of the coast, so first try to go right (to the W)
561 nXSeaward = nX - 1;
562 nYSeaward = nY;
563 nSeawardNewDirection = WEST;
564
565 // If can't do this, try to go straight on (to the S)
566 nXStraightOn = nX;
567 nYStraightOn = nY + 1;
568
569 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (E)
570 nXAntiSeaward = nX + 1;
571 nYAntiSeaward = nY;
572 nAntiSeawardNewDirection = EAST;
573
574 // As a last resort, go back (to the N)
575 nXGoBack = nX;
576 nYGoBack = nY - 1;
577 nGoBackNewDirection = NORTH;
578
579 break;
580
581 case WEST:
582 // The sea is towards the RHS (N) of the coast, so first try to go right (to the N)
583 nXSeaward = nX;
584 nYSeaward = nY - 1;
585 nSeawardNewDirection = NORTH;
586
587 // If can't do this, try to go straight on (to the W)
588 nXStraightOn = nX - 1;
589 nYStraightOn = nY;
590
591 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (S)
592 nXAntiSeaward = nX;
593 nYAntiSeaward = nY + 1;
594 nAntiSeawardNewDirection = SOUTH;
595
596 // As a last resort, go back (to the E)
597 nXGoBack = nX + 1;
598 nYGoBack = nY;
599 nGoBackNewDirection = EAST;
600
601 break;
602 }
603 break;
604
605 case LEFT_HANDED:
606 // The sea is to the left-hand side of the coast as we traverse it. We are just inland, so we need to keep heading left to find the sea
607 switch (nSearchDirection)
608 {
609 case NORTH:
610 // The sea is towards the LHS (W) of the coast, so first try to go left (to the W)
611 nXSeaward = nX - 1;
612 nYSeaward = nY;
613 nSeawardNewDirection = WEST;
614
615 // If can't do this, try to go straight on (to the N)
616 nXStraightOn = nX;
617 nYStraightOn = nY - 1;
618
619 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (E)
620 nXAntiSeaward = nX + 1;
621 nYAntiSeaward = nY;
622 nAntiSeawardNewDirection = EAST;
623
624 // As a last resort, go back (to the S)
625 nXGoBack = nX;
626 nYGoBack = nY + 1;
627 nGoBackNewDirection = SOUTH;
628
629 break;
630
631 case EAST:
632 // The sea is towards the LHS (N) of the coast, so first try to go left (to the N)
633 nXSeaward = nX;
634 nYSeaward = nY - 1;
635 nSeawardNewDirection = NORTH;
636
637 // If can't do this, try to go straight on (to the E)
638 nXStraightOn = nX + 1;
639 nYStraightOn = nY;
640
641 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (S)
642 nXAntiSeaward = nX;
643 nYAntiSeaward = nY + 1;
644 nAntiSeawardNewDirection = SOUTH;
645
646 // As a last resort, go back (to the W)
647 nXGoBack = nX - 1;
648 nYGoBack = nY;
649 nGoBackNewDirection = WEST;
650
651 break;
652
653 case SOUTH:
654 // The sea is towards the LHS (E) of the coast, so first try to go left (to the E)
655 nXSeaward = nX + 1;
656 nYSeaward = nY;
657 nSeawardNewDirection = EAST;
658
659 // If can't do this, try to go straight on (to the S)
660 nXStraightOn = nX;
661 nYStraightOn = nY + 1;
662
663 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (W)
664 nXAntiSeaward = nX - 1;
665 nYAntiSeaward = nY;
666 nAntiSeawardNewDirection = WEST;
667
668 // As a last resort, go back (to the N)
669 nXGoBack = nX;
670 nYGoBack = nY - 1;
671 nGoBackNewDirection = NORTH;
672
673 break;
674
675 case WEST:
676 // The sea is towards the LHS (S) of the coast, so first try to go left (to the S)
677 nXSeaward = nX;
678 nYSeaward = nY + 1;
679 nSeawardNewDirection = SOUTH;
680
681 // If can't do this, try to go straight on (to the W)
682 nXStraightOn = nX - 1;
683 nYStraightOn = nY;
684
685 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (N)
686 nXAntiSeaward = nX;
687 nYAntiSeaward = nY - 1;
688 nAntiSeawardNewDirection = NORTH;
689
690 // As a last resort, go back (to the E)
691 nXGoBack = nX + 1;
692 nYGoBack = nY;
693 nGoBackNewDirection = EAST;
694
695 break;
696 }
697 break;
698 }
699
700 // Now do the actual search for this timestep: first try going in the direction of the sea. Is this seaward cell still within the grid?
701 if (bIsWithinValidGrid(nXSeaward, nYSeaward))
702 {
703 // It is, so check if the cell in the seaward direction is a sea cell
704 if (m_pRasterGrid->m_Cell[nXSeaward][nYSeaward].bIsInContiguousSea())
705 {
706 // There is sea in this seaward direction, so we are on the coast
707 bAtCoast = true;
708
709 // Has the current cell already marked been marked as a coast cell?
710 if (! m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
711 {
712 // Not already marked, is this an intervention cell with the top above SWL?
713 if ((bIsInterventionCell(nX, nY)) && (m_pRasterGrid->m_Cell[nX][nY].dGetInterventionTopElev() >= m_dThisIterSWL))
714 {
715 // It is, so mark as coast and add it to the vector object
716 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
717 ILTempGridCRS.Append(&Pti);
718 }
719 else if (m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() >= m_dThisIterSWL)
720 {
721 // The sediment top is above SWL so mark as coast and add it to the vector object
722 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
723 ILTempGridCRS.Append(&Pti);
724 }
725 }
726 }
727 else
728 {
729 // The seaward cell is not a sea cell, so we will move to it next time
730 nX = nXSeaward;
731 nY = nYSeaward;
732
733 // And set a new search direction, to keep turning seaward
734 nSearchDirection = nSeawardNewDirection;
735 continue;
736 }
737 }
738
739 // OK, we couldn't move seaward (but we may have marked the current cell as coast) so next try to move straight on. Is this straight-ahead cell still within the grid?
740 if (bIsWithinValidGrid(nXStraightOn, nYStraightOn))
741 {
742 // It is, so check if there is sea immediately in front
743 if (m_pRasterGrid->m_Cell[nXStraightOn][nYStraightOn].bIsInContiguousSea())
744 {
745 // Sea is in front, so we are on the coast
746 bAtCoast = true;
747
748 // Has the current cell already marked been marked as a coast cell?
749 if (! m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
750 {
751 // Not already marked, is this an intervention cell with the top above SWL?
752 if ((bIsInterventionCell(nX, nY)) && (m_pRasterGrid->m_Cell[nX][nY].dGetInterventionTopElev() >= m_dThisIterSWL))
753 {
754 // It is, so mark as coast and add it to the vector object
755 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
756 ILTempGridCRS.Append(&Pti);
757 }
758 else if (m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() >= m_dThisIterSWL)
759 {
760 // The sediment top is above SWL so mark as coast and add it to the vector object
761 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
762 ILTempGridCRS.Append(&Pti);
763 }
764 }
765 }
766 else
767 {
768 // The straight-ahead cell is not a sea cell, so we will move to it next time
769 nX = nXStraightOn;
770 nY = nYStraightOn;
771
772 // The search direction remains unchanged
773 continue;
774 }
775 }
776
777 // Couldn't move either seaward or straight on (but we may have marked the current cell as coast) so next try to move in the anti-seaward direction. Is this anti-seaward cell still within the grid?
778 if (bIsWithinValidGrid(nXAntiSeaward, nYAntiSeaward))
779 {
780 // It is, so check if there is sea in this anti-seaward cell
781 if (m_pRasterGrid->m_Cell[nXAntiSeaward][nYAntiSeaward].bIsInContiguousSea())
782 {
783 // There is sea on the anti-seaward side, so we are on the coast
784 bAtCoast = true;
785
786 // Has the current cell already marked been marked as a coast cell?
787 if (! m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
788 {
789 // Not already marked, is this an intervention cell with the top above SWL?
790 if ((bIsInterventionCell(nX, nY)) && (m_pRasterGrid->m_Cell[nX][nY].dGetInterventionTopElev() >= m_dThisIterSWL))
791 {
792 // It is, so mark as coast and add it to the vector object
793 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
794 ILTempGridCRS.Append(&Pti);
795 }
796 else if (m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() >= m_dThisIterSWL)
797 {
798 // The sediment top is above SWL so mark as coast and add it to the vector object
799 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
800 ILTempGridCRS.Append(&Pti);
801 }
802 }
803 }
804 else
805 {
806 // The anti-seaward cell is not a sea cell, so we will move to it next time
807 nX = nXAntiSeaward;
808 nY = nYAntiSeaward;
809
810 // And set a new search direction, to keep turning seaward
811 nSearchDirection = nAntiSeawardNewDirection;
812 continue;
813 }
814 }
815
816 // Could not move to the seaward side, move straight ahead, or move to the anti-seaward side, so we must be in a single-cell dead end! As a last resort, turn round and move back to where we just came from, but first check that this is a valid cell
817 if (bIsWithinValidGrid(nXGoBack, nYGoBack))
818 {
819 nX = nXGoBack;
820 nY = nYGoBack;
821
822 // And change the search direction
823 nSearchDirection = nGoBackNewDirection;
824 }
825 else
826 {
827 // Our final choice is not a valid cell, so give up
828 bOffEdge = true;
829 break;
830 }
831 } while (true);
832
833 // OK, we have finished tracing this coastline on the grid. But is the coastline too long or too short?
834 int nCoastSize = ILTempGridCRS.nGetSize();
835
836 if (bOffEdge)
837 {
838 if (m_nLogFileDetail >= LOG_FILE_ALL)
839 LogStream << m_ulIter << ": ignoring temporary coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since hit off-edge cell at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}, coastline size is " << nCoastSize << endl;
840
841 // Unmark these cells as coast cells
842 for (int n = 0; n < nCoastSize; n++)
843 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
844
845 return RTN_ERR_TRACING_COAST;
846 }
847
848 if (bTooLong)
849 {
850 // Around loop too many times, so abandon this coastline
851 if (m_nLogFileDetail >= LOG_FILE_ALL)
852 {
853 LogStream << m_ulIter << ": ignoring temporary coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since round loop " << nRoundLoop << " times (m_nCoastMax = " << m_nCoastMax << "), coastline size is " << nCoastSize;
854
855 if (nCoastSize > 0)
856 LogStream << ", it ended at [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "}";
857
858 LogStream << endl;
859 }
860
861 // Unmark these cells as coast cells
862 for (int n = 0; n < nCoastSize; n++)
863 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
864
865 return RTN_ERR_TRACING_COAST;
866 }
867
868 if (bRepeating)
869 {
870 if (m_nLogFileDetail >= LOG_FILE_ALL)
871 {
872 LogStream << m_ulIter << ": ignoring temporary coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since repeating, coastline size is " << nCoastSize;
873
874 if (nCoastSize > 0)
875 LogStream << ", it ended at [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "}";
876
877 LogStream << endl;
878 }
879
880 // Unmark these cells as coast cells
881 for (int n = 0; n < nCoastSize; n++)
882 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
883
884 return RTN_ERR_TRACING_COAST;
885 }
886
887 if (nCoastSize == 0)
888 {
889 // Zero-length coastline, so abandon it
890 if (m_nLogFileDetail >= LOG_FILE_ALL)
891 LogStream << m_ulIter << ": abandoning zero-length coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "}" << endl;
892
893 return RTN_ERR_TRACING_COAST;
894 }
895
896 if (nCoastSize < m_nCoastMin)
897 {
898 // The vector coastline is too small, so abandon it
899 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
900 LogStream << m_ulIter << ": ignoring temporary coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} to [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "} since size (" << nCoastSize << ") is less than minimum (" << m_nCoastMin << ")" << endl;
901
902 // Unmark these cells as coast cells
903 for (int n = 0; n < nCoastSize; n++)
904 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
905
906 return RTN_ERR_TRACING_COAST;
907 }
908
909 // OK this new coastline is fine
910 int nEndX = nX;
911 int nEndY = nY;
912 int nCoastEndX = ILTempGridCRS[nCoastSize - 1].nGetX();
913 int nCoastEndY = ILTempGridCRS[nCoastSize - 1].nGetY();
914
915 if ((nCoastEndX != nEndX) || (nCoastEndY != nEndY))
916 {
917 // The grid-edge cell at nEndX, nEndY is not already at end of ILTempGridCRS. But is the final cell in ILTempGridCRS already at the edge of the grid?
918 if (! m_pRasterGrid->m_Cell[nCoastEndX][nCoastEndY].bIsBoundingBoxEdge())
919 {
920 // The final cell in ILTempGridCRS is not a grid-edge cell, so add the grid-edge cell and mark the cell as coastline
921 ILTempGridCRS.Append(nEndX, nEndY);
922 nCoastSize++;
923
924 m_pRasterGrid->m_Cell[nEndX][nEndY].SetAsCoastline(true);
925 }
926 }
927
928 // Need to specify start edge and end edge for smoothing routines
929 int nStartEdge = m_pRasterGrid->m_Cell[nStartX][nStartY].nGetBoundingBoxEdge();
930 int nEndEdge = m_pRasterGrid->m_Cell[nEndX][nEndY].nGetBoundingBoxEdge();
931
932 // Next, convert the grid coordinates in ILTempGridCRS (integer values stored as doubles) to external CRS coordinates (which will probably be non-integer, again stored as doubles). This is done now, so that smoothing is more effective
933 CGeomLine LTempExtCRS;
934 for (int j = 0; j < nCoastSize; j++)
935 {
936 LTempExtCRS.Append(dGridCentroidXToExtCRSX(ILTempGridCRS[j].nGetX()), dGridCentroidYToExtCRSY(ILTempGridCRS[j].nGetY()));
937 }
938
939 // Now do some smoothing of the vector output, if desired
940 if (m_nCoastSmooth == SMOOTH_RUNNING_MEAN)
941 LTempExtCRS = LSmoothCoastRunningMean(&LTempExtCRS);
942 else if (m_nCoastSmooth == SMOOTH_SAVITZKY_GOLAY)
943 LTempExtCRS = LSmoothCoastSavitzkyGolay(&LTempExtCRS, nStartEdge, nEndEdge);
944
945 // // DEBUG CODE ==================================================================================================
946 // LogStream << "==================================" << endl;
947 // for (int j = 0; j < nCoastSize; j++)
948 // {
949 // LogStream << "{" << dGridCentroidXToExtCRSX(ILTempGridCRS[j].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[j].nGetY()) << "}" << "\t{" << LTempExtCRS.dGetXAt(j) << ", " << LTempExtCRS.dGetYAt(j) << "}" << endl;
950 // }
951 // LogStream << "==================================" << endl;
952 // // DEBUG CODE ==================================================================================================
953
954 // Create a new coastline object and append to it the vector of coastline objects
955 CRWCoast CoastTmp(this);
956 m_VCoast.push_back(CoastTmp);
957 int nCoast = static_cast<int>(m_VCoast.size()) - 1;
958
959 // Set the coastline (Ext CRS)
960 m_VCoast[nCoast].SetCoastlineExtCRS(&LTempExtCRS);
961
962 // Set the coastline (Grid CRS)
963 m_VCoast[nCoast].SetCoastlineGridCRS(&ILTempGridCRS);
964
965 // CGeom2DPoint PtLast(DBL_MIN, DBL_MIN);
966 // for (int j = 0; j < nCoastSize; j++)
967 // {
968 // // Store the smoothed points (in external CRS) in the coast's m_LCoastlineExtCRS object, also append dummy values to the other attribute vectors
969 // if (PtLast != &LTempExtCRS[j]) // Avoid duplicate points
970 // {
971 // m_VCoast[nCoast].AppendPointToCoastlineExtCRS(LTempExtCRS[j].dGetX(), LTempExtCRS[j].dGetY());
972 //
973 // // Also store the locations of the corresponding unsmoothed points (in raster grid CRS) in the coast's m_ILCellsMarkedAsCoastline vector
974 // m_VCoast[nCoast].AppendCellMarkedAsCoastline(&ILTempGridCRS[j]);
975 // }
976 //
977 // PtLast = LTempExtCRS[j];
978 // }
979
980 // Set values for the coast's other attributes: set the coast's handedness, and start and end edges
981 m_VCoast[nCoast].SetSeaHandedness(nHandedness);
982 m_VCoast[nCoast].SetStartEdge(nStartEdge);
983 m_VCoast[nCoast].SetEndEdge(nEndEdge);
984
985 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
986 {
987 LogStream << m_ulIter << ": Coast " << nCoast << " created, from [" << nStartX << "][" << nStartY << "] to [" << nEndX << "][" << nEndY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} to {" << dGridCentroidXToExtCRSX(nEndX) << ", " << dGridCentroidYToExtCRSY(nEndY) << "} with " << nCoastSize << " points, handedness = " << (nHandedness == LEFT_HANDED ? "left" : "right") << endl;
988
989 LogStream << m_ulIter << ": Smoothed coastline " << nCoast << " runs from {" << LTempExtCRS[0].dGetX() << ", " << LTempExtCRS[0].dGetY() << "} to {" << LTempExtCRS[nCoastSize - 1].dGetX() << ", " << LTempExtCRS[nCoastSize - 1].dGetY() << "} i.e. from the ";
990 if (nStartEdge == NORTH)
991 LogStream << "north";
992 else if (nStartEdge == SOUTH)
993 LogStream << "south";
994 else if (nStartEdge == WEST)
995 LogStream << "west";
996 else if (nStartEdge == EAST)
997 LogStream << "east";
998 LogStream << " edge to the ";
999 if (nEndEdge == NORTH)
1000 LogStream << "north";
1001 else if (nEndEdge == SOUTH)
1002 LogStream << "south";
1003 else if (nEndEdge == WEST)
1004 LogStream << "west";
1005 else if (nEndEdge == EAST)
1006 LogStream << "east";
1007 LogStream << " edge" << endl;
1008 }
1009 // LogStream << "-----------------" << endl;
1010 // for (int kk = 0; kk < m_VCoast.back().nGetCoastlineSize(); kk++)
1011 // LogStream << kk << " [" << m_VCoast.back().pPtiGetCellMarkedAsCoastline(kk)->nGetX() << "][" << m_VCoast.back().pPtiGetCellMarkedAsCoastline(kk)->nGetY() << "] = {" << dGridCentroidXToExtCRSX(m_VCoast.back().pPtiGetCellMarkedAsCoastline(kk)->nGetX()) << ", " << dGridCentroidYToExtCRSY(m_VCoast.back().pPtiGetCellMarkedAsCoastline(kk)->nGetY()) << "}" << endl;
1012 // LogStream << "-----------------" << endl;
1013
1014 // Next calculate the curvature of the vector coastline
1015 DoCoastCurvature(nCoast, nHandedness);
1016
1017 // Calculate values for the coast's flux orientation vector
1018 CalcCoastTangents(nCoast);
1019
1020 // And create the vector of pointers to coastline-normal objects
1021 m_VCoast[nCoast].CreateProfilesAtCoastPoints();
1022
1023 return RTN_OK;
1024}
1025
1026//===============================================================================================================================
1028//===============================================================================================================================
1029int CSimulation::nLocateFloodAndCoasts(void)
1030{
1031 // Find all connected sea cells
1032 FindAllInundatedCells();
1033
1034 // Find every coastline on the raster grid, mark raster cells, then create the vector coastline
1035 int nRet = nTraceAllFloodCoasts();
1036 if (nRet != RTN_OK)
1037 return nRet;
1038
1039 // Have we created any coasts?
1040 switch (m_nLevel)
1041 {
1042 case 0: // WAVESETUP + SURGE:
1043 {
1044 if (m_VFloodWaveSetupSurge.empty())
1045 {
1046 cerr << m_ulIter << ": " << ERR << "no flood coastline located: this iteration SWL = " << m_dThisIterSWL << ", maximum DEM top surface elevation = " << m_dThisIterTopElevMax << ", minimum DEM top surface elevation = " << m_dThisIterTopElevMin << endl;
1047 return RTN_ERR_NOCOAST;
1048 }
1049 break;
1050 }
1051
1052 case 1: // WAVESETUP + SURGE + RUNUP:
1053 {
1054 if (m_VFloodWaveSetupSurgeRunup.empty())
1055 {
1056 cerr << m_ulIter << ": " << ERR << "no flood coastline located: this iteration SWL = " << m_dThisIterSWL << ", maximum DEM top surface elevation = " << m_dThisIterTopElevMax << ", minimum DEM top surface elevation = " << m_dThisIterTopElevMin << endl;
1057 return RTN_ERR_NOCOAST;
1058 }
1059 break;
1060 }
1061 }
1062
1063 return RTN_OK;
1064}
1065
1066//===============================================================================================================================
1068//===============================================================================================================================
1069int CSimulation::FindAllInundatedCells(void)
1070{
1071 for (int nX = 0; nX < m_nXGridSize; nX++)
1072 {
1073 for (int nY = 0; nY < m_nYGridSize; nY++)
1074 {
1075 m_pRasterGrid->m_Cell[nX][nY].UnSetCheckFloodCell(); // TODO 007 Do we need this?
1076 m_pRasterGrid->m_Cell[nX][nY].UnSetInContiguousFlood(); // TODO 007 Do we need this?
1077 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(false); // TODO 007 Do we need this?
1078 }
1079 }
1080 // Go along the list of edge cells
1081 for (unsigned int n = 0; n < m_VEdgeCell.size(); n++)
1082 {
1083 if (m_bOmitSearchNorthEdge && m_VEdgeCellEdge[n] == NORTH)
1084 continue;
1085
1086 if (m_bOmitSearchSouthEdge && m_VEdgeCellEdge[n] == SOUTH)
1087 continue;
1088
1089 if (m_bOmitSearchWestEdge && m_VEdgeCellEdge[n] == WEST)
1090 continue;
1091
1092 if (m_bOmitSearchEastEdge && m_VEdgeCellEdge[n] == EAST)
1093 continue;
1094
1095 int nX = m_VEdgeCell[n].nGetX();
1096 int nY = m_VEdgeCell[n].nGetY();
1097
1098 if ((! m_pRasterGrid->m_Cell[nX][nY].bIsCellFloodCheck()) && (m_pRasterGrid->m_Cell[nX][nY].bIsInundated()))
1099 {
1100 // This edge cell is below SWL and sea depth remains set to zero
1101 FloodFillLand(nX, nY);
1102 }
1103 }
1104
1105 return RTN_OK;
1106}
1107
1108//===============================================================================================================================
1110//===============================================================================================================================
1111void CSimulation::FloodFillLand(int const nXStart, int const nYStart)
1112{
1113 // The flood is at a user-specified location. So get the location from values read from the shapefile
1114 long unsigned int nLocIDs = m_VdFloodLocationX.size();
1115 double dDiffTotWaterLevel = 0;
1116
1117 double dAuxWaterLevelDiff = 0;
1118 if (nLocIDs == 0)
1119 {
1120 int pointCounter = 0;
1121 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
1122 {
1123 int nCoastSize = m_VCoast[nCoast].pLGetCoastlineExtCRS()->nGetSize();
1124 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
1125 {
1126 dAuxWaterLevelDiff = m_VCoast[nCoast].dGetLevel(nCoastPoint, m_nLevel);
1127 if (! isnan(dAuxWaterLevelDiff))
1128 {
1129 if (tAbs(dAuxWaterLevelDiff) < 1) // Limiting the maximum value that can be found (dAuxWaterLevelDiff != DBL_NODATA)
1130 {
1131 pointCounter++;
1132 dDiffTotWaterLevel += dAuxWaterLevelDiff;
1133 }
1134 }
1135 }
1136 }
1137 if (pointCounter > 0)
1138 dDiffTotWaterLevel /= pointCounter;
1139 else
1140 dDiffTotWaterLevel = 0;
1141 }
1142 else
1143 {
1144 for (long unsigned int n = 0; n < nLocIDs; n++)
1145 {
1146 double dPointGridXExtCRS = m_VdFloodLocationX[n];
1147 double dPointGridYExtCRS = m_VdFloodLocationY[n];
1148 double dMinDiffTotWaterLevelAtCoast = 1e10;
1149 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
1150 {
1151 int nCoastSize = m_VCoast[nCoast].pLGetCoastlineExtCRS()->nGetSize();
1152 double dMinDistSquare = 1e10;
1153 for (int nCoastPoint = 0; nCoastPoint < nCoastSize; nCoastPoint++)
1154 {
1155 double dCoastPointXExtCRS = m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint)->dGetX();
1156 double dCoastPointYExtCRS = m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(nCoastPoint)->dGetY();
1157
1158 double dDistSquare = (dCoastPointXExtCRS - dPointGridXExtCRS) * (dCoastPointXExtCRS - dPointGridXExtCRS) + (dCoastPointYExtCRS - dPointGridYExtCRS) * (dCoastPointYExtCRS - dPointGridYExtCRS);
1159
1160 if (dDistSquare < dMinDistSquare)
1161 {
1162 dAuxWaterLevelDiff = m_VCoast[nCoast].dGetLevel(nCoastPoint, m_nLevel);
1163 if (! isnan(dAuxWaterLevelDiff))
1164 {
1165 dMinDistSquare = dDistSquare;
1166 dMinDiffTotWaterLevelAtCoast = dAuxWaterLevelDiff;
1167 }
1168 }
1169 }
1170 }
1171 dDiffTotWaterLevel += dMinDiffTotWaterLevelAtCoast;
1172 }
1173 dDiffTotWaterLevel /= static_cast<double>(nLocIDs);
1174 }
1175
1176 m_dThisIterDiffTotWaterLevel = dDiffTotWaterLevel;
1177 switch (m_nLevel)
1178 {
1179 case 0: // WAVESETUP + STORMSURGE:
1180 m_dThisIterDiffWaveSetupSurgeWaterLevel = m_dThisIterDiffTotWaterLevel;
1181 break;
1182 case 1: // WAVESETUP + STORMSURGE + RUNUP:
1183 m_dThisIterDiffWaveSetupSurgeRunupWaterLevel = m_dThisIterDiffTotWaterLevel;
1184 break;
1185 }
1186
1187 // Create an empty stack
1188 stack<CGeom2DIPoint> PtiStackFlood;
1189
1190 // Start at the given edge cell, push this onto the stack
1191 PtiStackFlood.push(CGeom2DIPoint(nXStart, nYStart));
1192
1193 // Then do the flood fill loop until there are no more cell coordinates on the stack
1194 while (! PtiStackFlood.empty())
1195 {
1196 CGeom2DIPoint Pti = PtiStackFlood.top();
1197 PtiStackFlood.pop();
1198
1199 int nX = Pti.nGetX();
1200 int nY = Pti.nGetY();
1201
1202 while (nX >= 0)
1203 {
1204 if (m_pRasterGrid->m_Cell[nX][nY].bIsCellFloodCheck())
1205 break;
1206 if (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel())
1207 break;
1208 nX--;
1209 }
1210
1211 nX++;
1212
1213 bool bSpanAbove = false;
1214 bool bSpanBelow = false;
1215
1216 while (nX < m_nXGridSize)
1217 {
1218 if (m_pRasterGrid->m_Cell[nX][nY].bIsCellFloodCheck())
1219 break;
1220 if (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel())
1221 break;
1222
1223 // Flood this cell
1224 m_pRasterGrid->m_Cell[nX][nY].SetCheckFloodCell(); // TODO 007 Do we need this?
1225 m_pRasterGrid->m_Cell[nX][nY].SetInContiguousFlood(); // TODO 007 Do we need this?
1226
1227 switch (m_nLevel)
1228 {
1229 case 0: // WAVESETUP + STORMSURGE:
1230 m_pRasterGrid->m_Cell[nX][nY].SetFloodBySetupSurge();
1231 break;
1232 case 1: // WAVESETUP + STORMSURGE + RUNUP:
1233 m_pRasterGrid->m_Cell[nX][nY].SetFloodBySetupSurgeRunup();
1234 break;
1235 }
1236
1237 if ((! bSpanAbove) && (nY > 0) && (m_pRasterGrid->m_Cell[nX][nY - 1].bIsElevLessThanWaterLevel()) && (! m_pRasterGrid->m_Cell[nX][nY - 1].bIsCellFloodCheck()))
1238 {
1239 PtiStackFlood.push(CGeom2DIPoint(nX, nY - 1));
1240 bSpanAbove = true;
1241 }
1242 else if (bSpanAbove && (nY > 0) && (! m_pRasterGrid->m_Cell[nX][nY - 1].bIsElevLessThanWaterLevel()))
1243 {
1244 bSpanAbove = false;
1245 }
1246
1247 if ((! bSpanBelow) && (nY < m_nYGridSize - 1) && (m_pRasterGrid->m_Cell[nX][nY + 1].bIsElevLessThanWaterLevel()) && (! m_pRasterGrid->m_Cell[nX][nY + 1].bIsCellFloodCheck()))
1248 {
1249 PtiStackFlood.push(CGeom2DIPoint(nX, nY + 1));
1250 bSpanBelow = true;
1251 }
1252 else if (bSpanBelow && (nY < m_nYGridSize - 1) && (! m_pRasterGrid->m_Cell[nX][nY + 1].bIsElevLessThanWaterLevel()))
1253 {
1254 bSpanBelow = false;
1255 }
1256
1257 nX++;
1258 }
1259 }
1260}
1261
1262//===============================================================================================================================
1264//===============================================================================================================================
1265int CSimulation::nTraceAllFloodCoasts(void)
1266{
1267 vector<bool> VbPossibleStartCellLHEdge;
1268 vector<bool> VbTraced;
1269 vector<int> VnSearchDirection;
1270 vector<CGeom2DIPoint> V2DIPossibleStartCell;
1271
1272 // Go along the list of edge cells and look for possible coastline start cells
1273 for (unsigned int n = 0; n < m_VEdgeCell.size() - 1; n++)
1274 {
1275 if (m_bOmitSearchNorthEdge && (m_VEdgeCellEdge[n] == NORTH || m_VEdgeCellEdge[n + 1] == NORTH))
1276 continue;
1277
1278 if (m_bOmitSearchSouthEdge && (m_VEdgeCellEdge[n] == SOUTH || m_VEdgeCellEdge[n + 1] == SOUTH))
1279 continue;
1280
1281 if (m_bOmitSearchWestEdge && (m_VEdgeCellEdge[n] == WEST || m_VEdgeCellEdge[n + 1] == WEST))
1282 continue;
1283
1284 if (m_bOmitSearchEastEdge && (m_VEdgeCellEdge[n] == EAST || m_VEdgeCellEdge[n + 1] == EAST))
1285 continue;
1286
1287 int nXThis = m_VEdgeCell[n].nGetX();
1288 int nYThis = m_VEdgeCell[n].nGetY();
1289 int nXNext = m_VEdgeCell[n + 1].nGetX();
1290 int nYNext = m_VEdgeCell[n + 1].nGetY();
1291
1292 // Get "Is it sea?" information for 'this' and 'next' cells TODO 007 Not clear
1293 bool bThisCellIsSea = m_pRasterGrid->m_Cell[nXThis][nYThis].bIsInContiguousFlood(); // TODO 007 Do we need this?
1294 bool bNextCellIsSea = m_pRasterGrid->m_Cell[nXNext][nYNext].bIsInContiguousFlood(); // TODO 007 Do we need this?
1295
1296 // Are we at a coast?
1297 if ((! bThisCellIsSea) && bNextCellIsSea)
1298 {
1299 // 'This' cell is just inland, has it already been flagged as a possible start for a coastline (even if this subsequently 'failed' as a coastline)?
1300 // if (! m_pRasterGrid->m_Cell[nXThis][nYThis].bIsPossibleCoastStartCell())
1301 {
1302 // It has not, so flag it
1303 m_pRasterGrid->m_Cell[nXThis][nYThis].SetPossibleFloodStartCell();
1304
1305 // And save it
1306 V2DIPossibleStartCell.push_back(CGeom2DIPoint(nXThis, nYThis));
1307 VbPossibleStartCellLHEdge.push_back(true);
1308 VnSearchDirection.push_back(nGetOppositeDirection(m_VEdgeCellEdge[n]));
1309 VbTraced.push_back(false);
1310 }
1311 }
1312 else if (bThisCellIsSea && (! bNextCellIsSea))
1313 {
1314 // The 'next' cell is just inland, has it already been flagged as a possible start for a coastline (even if this subsequently 'failed' as a coastline)?
1315 // if (! m_pRasterGrid->m_Cell[nXNext][nYNext].bIsPossibleCoastStartCell())
1316 {
1317 // It has not, so flag it
1318 m_pRasterGrid->m_Cell[nXNext][nYNext].SetPossibleFloodStartCell();
1319
1320 // And save it
1321 V2DIPossibleStartCell.push_back(CGeom2DIPoint(nXNext, nYNext));
1322 VbPossibleStartCellLHEdge.push_back(false);
1323 VnSearchDirection.push_back(nGetOppositeDirection(m_VEdgeCellEdge[n + 1]));
1324 VbTraced.push_back(false);
1325 }
1326 }
1327 }
1328
1329 bool bAtLeastOneCoastTraced = false;
1330 for (unsigned int n = 0; n < V2DIPossibleStartCell.size(); n++)
1331 {
1332 if (! VbTraced[n])
1333 {
1334 int nRet = 0;
1335 if (VbPossibleStartCellLHEdge[n])
1336 {
1337 nRet = nTraceFloodCoastLine(n, VnSearchDirection[n], LEFT_HANDED, &VbTraced, &V2DIPossibleStartCell);
1338 }
1339 else
1340 {
1341 nRet = nTraceFloodCoastLine(n, VnSearchDirection[n], RIGHT_HANDED, &VbTraced, &V2DIPossibleStartCell);
1342 }
1343
1344 if (nRet == RTN_OK)
1345 {
1346 // We have a valid coastline starting from this possible start cell
1347 VbTraced[n] = true;
1348 bAtLeastOneCoastTraced = true;
1349 }
1350 }
1351 }
1352
1353 if (bAtLeastOneCoastTraced)
1354 return RTN_OK;
1355 else
1356 return RTN_ERR_TRACING_COAST;
1357}
1358
1359//===============================================================================================================================
1361//===============================================================================================================================
1362int CSimulation::nTraceFloodCoastLine(unsigned int const nTraceFromStartCellIndex, int const nStartSearchDirection, int const nHandedness, vector<bool>* pVbTraced, vector<CGeom2DIPoint> const* pV2DIPossibleStartCell)
1363{
1364 bool bHitStartCell = false;
1365 bool bAtCoast = false;
1366 bool bHasLeftStartEdge = false;
1367 bool bTooLong = false;
1368 bool bOffEdge = false;
1369 bool bRepeating = false;
1370
1371 int nStartX = pV2DIPossibleStartCell->at(nTraceFromStartCellIndex).nGetX();
1372 int nStartY = pV2DIPossibleStartCell->at(nTraceFromStartCellIndex).nGetY();
1373 int nX = nStartX;
1374 int nY = nStartY;
1375 int nSearchDirection = nStartSearchDirection;
1376 int nRoundLoop = -1;
1377 // nThisLen = 0;
1378 // nLastLen = 0,
1379 // nPreLastLen = 0;
1380
1381 // Temporary coastline as integer points (grid CRS)
1382 CGeomILine ILTempGridCRS;
1383
1384 // Mark the start cell as coast and add it to the vector object
1385 m_pRasterGrid->m_Cell[nStartX][nStartY].SetAsFloodLine(true);
1386 CGeom2DIPoint PtiStart(nStartX, nStartY);
1387 ILTempGridCRS.Append(&PtiStart);
1388
1389 // Start at this grid-edge point and trace the rest of the coastline using the 'wall follower' rule for maze traversal, trying to keep next to cells flagged as sea
1390 do
1391 {
1392 // // DEBUG CODE ==============================================================================================================
1393 // LogStream << "Now at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
1394 // LogStream << "ILTempGridCRS is now:" << endl;
1395 // for (int n = 0; n < ILTempGridCRS.nGetSize(); n++)
1396 // LogStream << "[" << ILTempGridCRS[n].nGetX() << "][" << ILTempGridCRS[n].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[n].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[n].nGetY()) << "}" << endl;
1397 // LogStream << "=================" << endl;
1398 // // DEBUG CODE ==============================================================================================================
1399
1400 // Safety check
1401 if (++nRoundLoop > m_nCoastMax)
1402 {
1403 bTooLong = true;
1404
1405 // LogStream << m_ulIter << ": abandoning coastline tracing from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "}, exceeded maximum search length (" << m_nCoastMax << ")" << endl;
1406
1407 // for (int n = 0; n < ILTempGridCRS.nGetSize(); n++)
1408 // LogStream << "[" << ILTempGridCRS[n].nGetX() << "][" << ILTempGridCRS[n].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[n].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[n].nGetY()) << "}" << endl;
1409 // LogStream << endl;
1410
1411 break;
1412 }
1413
1414 // Another safety check
1415 if ((nRoundLoop > 10) && (ILTempGridCRS.nGetSize() < 2))
1416 {
1417 // We've been 10 times round the loop but the coast is still less than 2 coastline points in length, so we must be repeating
1418 bRepeating = true;
1419
1420 break;
1421 }
1422
1423 // OK so far: so have we left the start edge?
1424 if (! bHasLeftStartEdge)
1425 {
1426 // We have not yet left the start edge
1427 if (((nStartSearchDirection == SOUTH) && (nY > nStartY)) || ((nStartSearchDirection == NORTH) && (nY < nStartY)) ||
1428 ((nStartSearchDirection == EAST) && (nX > nStartX)) || ((nStartSearchDirection == WEST) && (nX < nStartX)))
1429 bHasLeftStartEdge = true;
1430
1431 // Flag this cell to ensure that it is not chosen as a coastline start cell later
1432 m_pRasterGrid->m_Cell[nX][nY].SetPossibleFloodStartCell();
1433 // LogStream << "Flagging [" << nX << "][" << nY << "] as possible coast start cell NOT YET LEFT EDGE" << endl;
1434 }
1435
1436 // Leave the loop if the vector coastline has left the start edge, then we find a coast cell which is a possible start cell from which a coastline has not yet been traced
1437 // if (bHasLeftStartEdge && bAtCoast)
1438 {
1439 for (unsigned int nn = 0; nn < pVbTraced->size(); nn++)
1440 {
1441 if ((nn != nTraceFromStartCellIndex) && (!pVbTraced->at(nn)))
1442 {
1443 // LogStream << "[" << pV2DIPossibleStartCell->at(nn).nGetX() << "][" << pV2DIPossibleStartCell->at(nn).nGetY() << "]" << endl;
1444
1445 if (bAtCoast && (nX == pV2DIPossibleStartCell->at(nn).nGetX()) && (nY == pV2DIPossibleStartCell->at(nn).nGetY()))
1446 {
1447 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
1448 LogStream << m_ulIter << ": valid flood coastline found, traced from [" << nStartX << "][" << nStartY << "] and hit another start cell at [" << nX << "][" << nY << "]" << endl;
1449
1450 pVbTraced->at(nn) = true;
1451 bHitStartCell = true;
1452 break;
1453 }
1454 }
1455 }
1456 // LogStream << endl;
1457 }
1458
1459 if (bHitStartCell)
1460 break;
1461
1462 // OK now sort out the next iteration of the search
1463 int nXSeaward = 0;
1464 int nYSeaward = 0;
1465 int nSeawardNewDirection = 0;
1466 int nXStraightOn = 0;
1467 int nYStraightOn = 0;
1468 int nXAntiSeaward = 0;
1469 int nYAntiSeaward = 0;
1470 int nAntiSeawardNewDirection = 0;
1471 int nXGoBack = 0;
1472 int nYGoBack = 0;
1473 int nGoBackNewDirection = 0;
1474
1475 CGeom2DIPoint Pti(nX, nY);
1476
1477 // Set up the variables
1478 switch (nHandedness)
1479 {
1480 case RIGHT_HANDED:
1481 // The sea is to the right-hand side of the coast as we traverse it. We are just inland, so we need to keep heading right to find the sea
1482 switch (nSearchDirection)
1483 {
1484 case NORTH:
1485 // The sea is towards the RHS (E) of the coast, so first try to go right (to the E)
1486 nXSeaward = nX + 1;
1487 nYSeaward = nY;
1488 nSeawardNewDirection = EAST;
1489
1490 // If can't do this, try to go straight on (to the N)
1491 nXStraightOn = nX;
1492 nYStraightOn = nY - 1;
1493
1494 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (W)
1495 nXAntiSeaward = nX - 1;
1496 nYAntiSeaward = nY;
1497 nAntiSeawardNewDirection = WEST;
1498
1499 // As a last resort, go back (to the S)
1500 nXGoBack = nX;
1501 nYGoBack = nY + 1;
1502 nGoBackNewDirection = SOUTH;
1503
1504 break;
1505
1506 case EAST:
1507 // The sea is towards the RHS (S) of the coast, so first try to go right (to the S)
1508 nXSeaward = nX;
1509 nYSeaward = nY + 1;
1510 nSeawardNewDirection = SOUTH;
1511
1512 // If can't do this, try to go straight on (to the E)
1513 nXStraightOn = nX + 1;
1514 nYStraightOn = nY;
1515
1516 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (N)
1517 nXAntiSeaward = nX;
1518 nYAntiSeaward = nY - 1;
1519 nAntiSeawardNewDirection = NORTH;
1520
1521 // As a last resort, go back (to the W)
1522 nXGoBack = nX - 1;
1523 nYGoBack = nY;
1524 nGoBackNewDirection = WEST;
1525
1526 break;
1527
1528 case SOUTH:
1529 // The sea is towards the RHS (W) of the coast, so first try to go right (to the W)
1530 nXSeaward = nX - 1;
1531 nYSeaward = nY;
1532 nSeawardNewDirection = WEST;
1533
1534 // If can't do this, try to go straight on (to the S)
1535 nXStraightOn = nX;
1536 nYStraightOn = nY + 1;
1537
1538 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (E)
1539 nXAntiSeaward = nX + 1;
1540 nYAntiSeaward = nY;
1541 nAntiSeawardNewDirection = EAST;
1542
1543 // As a last resort, go back (to the N)
1544 nXGoBack = nX;
1545 nYGoBack = nY - 1;
1546 nGoBackNewDirection = NORTH;
1547
1548 break;
1549
1550 case WEST:
1551 // The sea is towards the RHS (N) of the coast, so first try to go right (to the N)
1552 nXSeaward = nX;
1553 nYSeaward = nY - 1;
1554 nSeawardNewDirection = NORTH;
1555
1556 // If can't do this, try to go straight on (to the W)
1557 nXStraightOn = nX - 1;
1558 nYStraightOn = nY;
1559
1560 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (S)
1561 nXAntiSeaward = nX;
1562 nYAntiSeaward = nY + 1;
1563 nAntiSeawardNewDirection = SOUTH;
1564
1565 // As a last resort, go back (to the E)
1566 nXGoBack = nX + 1;
1567 nYGoBack = nY;
1568 nGoBackNewDirection = EAST;
1569
1570 break;
1571 }
1572 break;
1573
1574 case LEFT_HANDED:
1575 // The sea is to the left-hand side of the coast as we traverse it. We are just inland, so we need to keep heading left to find the sea
1576 switch (nSearchDirection)
1577 {
1578 case NORTH:
1579 // The sea is towards the LHS (W) of the coast, so first try to go left (to the W)
1580 nXSeaward = nX - 1;
1581 nYSeaward = nY;
1582 nSeawardNewDirection = WEST;
1583
1584 // If can't do this, try to go straight on (to the N)
1585 nXStraightOn = nX;
1586 nYStraightOn = nY - 1;
1587
1588 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (E)
1589 nXAntiSeaward = nX + 1;
1590 nYAntiSeaward = nY;
1591 nAntiSeawardNewDirection = EAST;
1592
1593 // As a last resort, go back (to the S)
1594 nXGoBack = nX;
1595 nYGoBack = nY + 1;
1596 nGoBackNewDirection = SOUTH;
1597
1598 break;
1599
1600 case EAST:
1601 // The sea is towards the LHS (N) of the coast, so first try to go left (to the N)
1602 nXSeaward = nX;
1603 nYSeaward = nY - 1;
1604 nSeawardNewDirection = NORTH;
1605
1606 // If can't do this, try to go straight on (to the E)
1607 nXStraightOn = nX + 1;
1608 nYStraightOn = nY;
1609
1610 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (S)
1611 nXAntiSeaward = nX;
1612 nYAntiSeaward = nY + 1;
1613 nAntiSeawardNewDirection = SOUTH;
1614
1615 // As a last resort, go back (to the W)
1616 nXGoBack = nX - 1;
1617 nYGoBack = nY;
1618 nGoBackNewDirection = WEST;
1619
1620 break;
1621
1622 case SOUTH:
1623 // The sea is towards the LHS (E) of the coast, so first try to go left (to the E)
1624 nXSeaward = nX + 1;
1625 nYSeaward = nY;
1626 nSeawardNewDirection = EAST;
1627
1628 // If can't do this, try to go straight on (to the S)
1629 nXStraightOn = nX;
1630 nYStraightOn = nY + 1;
1631
1632 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (W)
1633 nXAntiSeaward = nX - 1;
1634 nYAntiSeaward = nY;
1635 nAntiSeawardNewDirection = WEST;
1636
1637 // As a last resort, go back (to the N)
1638 nXGoBack = nX;
1639 nYGoBack = nY - 1;
1640 nGoBackNewDirection = NORTH;
1641
1642 break;
1643
1644 case WEST:
1645 // The sea is towards the LHS (S) of the coast, so first try to go left (to the S)
1646 nXSeaward = nX;
1647 nYSeaward = nY + 1;
1648 nSeawardNewDirection = SOUTH;
1649
1650 // If can't do this, try to go straight on (to the W)
1651 nXStraightOn = nX - 1;
1652 nYStraightOn = nY;
1653
1654 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (N)
1655 nXAntiSeaward = nX;
1656 nYAntiSeaward = nY - 1;
1657 nAntiSeawardNewDirection = NORTH;
1658
1659 // As a last resort, go back (to the E)
1660 nXGoBack = nX + 1;
1661 nYGoBack = nY;
1662 nGoBackNewDirection = EAST;
1663
1664 break;
1665 }
1666 break;
1667 }
1668
1669 // Now do the actual search for this timestep: first try going in the direction of the sea. Is this seaward cell still within the grid?
1670 if (bIsWithinValidGrid(nXSeaward, nYSeaward))
1671 {
1672 // It is, so check if the cell in the seaward direction is a sea cell TODO 007 Not clear
1673 if (m_pRasterGrid->m_Cell[nXSeaward][nYSeaward].bIsInContiguousFlood()) // TODO 007 Do we need this?
1674 {
1675 // There is sea in this seaward direction, so we are on the coast
1676 bAtCoast = true;
1677
1678 // Has the current cell already marked been marked as a coast cell?
1679 if (! m_pRasterGrid->m_Cell[nX][nY].bIsFloodLine())
1680 {
1681 // Not already marked, is this an intervention cell with the top above SWL?
1682 if ((bIsInterventionCell(nX, nY)) && (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel()))
1683 {
1684 // It is, so mark as coast and add it to the vector object
1685 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(true);
1686 ILTempGridCRS.Append(&Pti);
1687 }
1688 else if (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel())
1689 {
1690 // The sediment top is above SWL so mark as coast and add it to the vector object
1691 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(true);
1692 ILTempGridCRS.Append(&Pti);
1693 }
1694 }
1695 }
1696 else
1697 {
1698 // The seaward cell is not a sea cell, so we will move to it next time
1699 nX = nXSeaward;
1700 nY = nYSeaward;
1701
1702 // And set a new search direction, to keep turning seaward
1703 nSearchDirection = nSeawardNewDirection;
1704 continue;
1705 }
1706 }
1707
1708 // OK, we couldn't move seaward (but we may have marked the current cell as coast) so next try to move straight on. Is this straight-ahead cell still within the grid?
1709 if (bIsWithinValidGrid(nXStraightOn, nYStraightOn))
1710 {
1711 // It is, so check if there is sea immediately in front TODO 007 Not clear
1712 if (m_pRasterGrid->m_Cell[nXStraightOn][nYStraightOn].bIsInContiguousFlood()) // TODO 007 Do we need this?
1713 {
1714 // Sea is in front, so we are on the coast
1715 bAtCoast = true;
1716
1717 // Has the current cell already marked been marked as a coast cell?
1718 if (! m_pRasterGrid->m_Cell[nX][nY].bIsFloodLine())
1719 {
1720 // Not already marked, is this an intervention cell with the top above SWL?
1721 if ((bIsInterventionCell(nX, nY)) && (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel()))
1722 {
1723 // It is, so mark as coast and add it to the vector object
1724 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(true);
1725 ILTempGridCRS.Append(&Pti);
1726 }
1727 else if (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel())
1728 {
1729 // The sediment top is above SWL so mark as coast and add it to the vector object
1730 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(true);
1731 ILTempGridCRS.Append(&Pti);
1732 }
1733 }
1734 }
1735 else
1736 {
1737 // The straight-ahead cell is not a sea cell, so we will move to it next time
1738 nX = nXStraightOn;
1739 nY = nYStraightOn;
1740
1741 // The search direction remains unchanged
1742 continue;
1743 }
1744 }
1745
1746 // Couldn't move either seaward or straight on (but we may have marked the current cell as coast) so next try to move in the anti-seaward direction. Is this anti-seaward cell still within the grid?
1747 if (bIsWithinValidGrid(nXAntiSeaward, nYAntiSeaward))
1748 {
1749 // It is, so check if there is sea in this anti-seaward cell TODO 007 Not clear
1750 if (m_pRasterGrid->m_Cell[nXAntiSeaward][nYAntiSeaward].bIsInContiguousFlood()) // TODO 007 Do we need this?
1751 {
1752 // There is sea on the anti-seaward side, so we are on the coast
1753 bAtCoast = true;
1754
1755 // Has the current cell already marked been marked as a coast cell?
1756 if (! m_pRasterGrid->m_Cell[nX][nY].bIsFloodLine())
1757 {
1758 // Not already marked, is this an intervention cell with the top above SWL?
1759 if ((bIsInterventionCell(nX, nY)) && (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel()))
1760 {
1761 // It is, so mark as coast and add it to the vector object
1762 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(true);
1763 ILTempGridCRS.Append(&Pti);
1764 }
1765 else if (! m_pRasterGrid->m_Cell[nX][nY].bIsElevLessThanWaterLevel())
1766 {
1767 // The sediment top is above SWL so mark as coast and add it to the vector object
1768 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(true);
1769 ILTempGridCRS.Append(&Pti);
1770 }
1771 }
1772 }
1773 else
1774 {
1775 // The anti-seaward cell is not a sea cell, so we will move to it next time
1776 nX = nXAntiSeaward;
1777 nY = nYAntiSeaward;
1778
1779 // And set a new search direction, to keep turning seaward
1780 nSearchDirection = nAntiSeawardNewDirection;
1781 continue;
1782 }
1783 }
1784
1785 // Could not move to the seaward side, move straight ahead, or move to the anti-seaward side, so we must be in a single-cell dead end! As a last resort, turn round and move back to where we just came from, but first check that this is a valid cell
1786 if (bIsWithinValidGrid(nXGoBack, nYGoBack))
1787 {
1788 nX = nXGoBack;
1789 nY = nYGoBack;
1790
1791 // And change the search direction
1792 nSearchDirection = nGoBackNewDirection;
1793 }
1794 else
1795 {
1796 // Our final choice is not a valid cell, so give up
1797 bOffEdge = true;
1798 break;
1799 }
1800 } while (true);
1801
1802 // OK, we have finished tracing this coastline on the grid. But is the coastline too long or too short?
1803 int nCoastSize = ILTempGridCRS.nGetSize();
1804
1805 if (bOffEdge)
1806 {
1807 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
1808 LogStream << m_ulIter << ": ignoring temporary flood coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since hit off-edge cell at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}, coastline size is " << nCoastSize << endl;
1809
1810 // Unmark these cells as coast cells
1811 for (int n = 0; n < nCoastSize; n++)
1812 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsFloodLine(false);
1813
1814 return RTN_ERR_TRACING_COAST;
1815 }
1816
1817 if (bTooLong)
1818 {
1819 // Around loop too many times, so abandon this coastline
1820 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
1821 {
1822 LogStream << m_ulIter << ": ignoring temporary flood coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since round loop " << nRoundLoop << " times (m_nCoastMax = " << m_nCoastMax << "), coastline size is " << nCoastSize;
1823
1824 if (nCoastSize > 0)
1825 LogStream << ", it ended at [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "}";
1826
1827 LogStream << endl;
1828 }
1829
1830 // Unmark these cells as coast cells
1831 for (int n = 0; n < nCoastSize; n++)
1832 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsFloodLine(false);
1833
1834 return RTN_ERR_TRACING_COAST;
1835 }
1836
1837 if (bRepeating)
1838 {
1839 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
1840 {
1841 LogStream << m_ulIter << ": ignoring temporary flood coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since repeating, coastline size is " << nCoastSize;
1842
1843 if (nCoastSize > 0)
1844 LogStream << ", it ended at [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "}";
1845
1846 LogStream << endl;
1847 }
1848
1849 // Unmark these cells as coast cells
1850 for (int n = 0; n < nCoastSize; n++)
1851 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsFloodLine(false);
1852
1853 return RTN_ERR_TRACING_COAST;
1854 }
1855
1856 if (nCoastSize == 0)
1857 {
1858 // Zero-length coastline, so abandon it
1859 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
1860 LogStream << m_ulIter << ": abandoning zero-length flood coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "}" << endl;
1861
1862 return RTN_ERR_TRACING_COAST;
1863 }
1864
1865 if (nCoastSize < m_nCoastMin)
1866 {
1867 // The vector coastline is too small, so abandon it
1868 if (m_nLogFileDetail >= LOG_FILE_HIGH_DETAIL)
1869 LogStream << m_ulIter << ": ignoring temporary flood coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} to [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "} since size (" << nCoastSize << ") is less than minimum (" << m_nCoastMin << ")" << endl;
1870
1871 // Unmark these cells as coast cells
1872 for (int n = 0; n < nCoastSize; n++)
1873 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsFloodLine(false);
1874
1875 return RTN_ERR_TRACING_COAST;
1876 }
1877
1878 // OK this new coastline is fine
1879 int nEndX = nX;
1880 int nEndY = nY;
1881 int nCoastEndX = ILTempGridCRS[nCoastSize - 1].nGetX();
1882 int nCoastEndY = ILTempGridCRS[nCoastSize - 1].nGetY();
1883
1884 if ((nCoastEndX != nEndX) || (nCoastEndY != nEndY))
1885 {
1886 // The grid-edge cell at nEndX, nEndY is not already at end of ILTempGridCRS. But is the final cell in ILTempGridCRS already at the edge of the grid?
1887 if (! m_pRasterGrid->m_Cell[nCoastEndX][nCoastEndY].bIsBoundingBoxEdge())
1888 {
1889 // The final cell in ILTempGridCRS is not a grid-edge cell, so add the grid-edge cell and mark the cell as coastline
1890 ILTempGridCRS.Append(nEndX, nEndY);
1891 nCoastSize++;
1892
1893 m_pRasterGrid->m_Cell[nEndX][nEndY].SetAsFloodLine(true);
1894 }
1895 }
1896
1897 // Need to specify start edge and end edge for smoothing routines
1898 // int
1899 // nStartEdge = m_pRasterGrid->m_Cell[nStartX][nStartY].nGetBoundingBoxEdge(),
1900 // nEndEdge = m_pRasterGrid->m_Cell[nEndX][nEndY].nGetBoundingBoxEdge();
1901
1902 // Next, convert the grid coordinates in ILTempGridCRS (integer values stored as doubles) to external CRS coordinates (which will probably be non-integer, again stored as doubles). This is done now, so that smoothing is more effective
1903 CGeomLine LTempExtCRS;
1904 for (int j = 0; j < nCoastSize; j++)
1905 {
1906 LTempExtCRS.Append(dGridCentroidXToExtCRSX(ILTempGridCRS[j].nGetX()), dGridCentroidYToExtCRSY(ILTempGridCRS[j].nGetY()));
1907 }
1908
1909 // Now do some smoothing of the vector output, if desired
1910 // if (m_nCoastSmooth == SMOOTH_RUNNING_MEAN)
1911 // LTempExtCRS = LSmoothCoastRunningMean(&LTempExtCRS);
1912 // else if (m_nCoastSmooth == SMOOTH_SAVITZKY_GOLAY)
1913 // LTempExtCRS = LSmoothCoastSavitzkyGolay(&LTempExtCRS, nStartEdge, nEndEdge);
1914
1915 // // DEBUG CODE ==================================================================================================
1916 // LogStream << "==================================" << endl;
1917 // for (int j = 0; j < nCoastSize; j++)
1918 // {
1919 // LogStream << "{" << dGridCentroidXToExtCRSX(ILTempGridCRS[j].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[j].nGetY()) << "}" << "\t{" << LTempExtCRS.dGetXAt(j) << ", " << LTempExtCRS.dGetYAt(j) << "}" << endl;
1920 // }
1921 // LogStream << "==================================" << endl;
1922 // // DEBUG CODE ==================================================================================================
1923
1924 // Create a new coastline object and append to it the vector of coastline objects
1925 CRWCoast CoastTmp(this);
1926 int nCoast;
1927 switch (m_nLevel)
1928 {
1929 case 0:
1930 m_VFloodWaveSetupSurge.push_back(CoastTmp);
1931 nCoast = static_cast<int>(m_VFloodWaveSetupSurge.size()) - 1;
1932 m_VFloodWaveSetupSurge[nCoast].SetCoastlineExtCRS(&LTempExtCRS);
1933 break;
1934
1935 case 1:
1936 m_VFloodWaveSetupSurgeRunup.push_back(CoastTmp);
1937 nCoast = static_cast<int>(m_VFloodWaveSetupSurgeRunup.size()) - 1;
1938 m_VFloodWaveSetupSurgeRunup[nCoast].SetCoastlineExtCRS(&LTempExtCRS);
1939 }
1940
1941 return RTN_OK;
1942}
CA2DIShape::nGetSize
int nGetSize(void) const
Returns the number of integer point in the vector which represents this 2D shape.
Definition 2di_shape.cpp:70
CA2DIShape::Append
void Append(CGeom2DIPoint const *)
Appends a new integer point to the vector which represents this 2D shape.
Definition 2di_shape.cpp:81
CA2DShape::Append
void Append(CGeom2DPoint const *)
Appends a point to this 2D shape.
Definition 2d_shape.cpp:67
CGeom2DIPoint
Geometry class used to represent 2D point objects with integer coordinates.
Definition 2di_point.h:29
CGeom2DIPoint::nGetY
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
Definition 2di_point.cpp:48
CGeom2DIPoint::nGetX
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Definition 2di_point.cpp:42
CGeomILine
Geometry class used to represent 2D vector integer line objects.
Definition i_line.h:31
CGeomLine
Geometry class used to represent 2D vector line objects.
Definition line.h:31
CRWCellLandform
Real-world class used to represent the landform of a cell.
Definition cell_landform.h:29
CRWCellLandform::nGetLFCategory
int nGetLFCategory(void) const
Get the landform category.
Definition cell_landform.cpp:49
CRWCellLandform::SetLFCategory
void SetLFCategory(int const)
Set the landform category.
Definition cell_landform.cpp:43
CRWCoast
Real-world class used to represent coastline objects.
Definition coast.h:39
CSimulation::m_nLogFileDetail
int m_nLogFileDetail
The level of detail in the log file output. Can be LOG_FILE_LOW_DETAIL, LOG_FILE_MIDDLE_DETAIL,...
Definition simulation.h:530
CSimulation::m_nYMinBoundingBox
int m_nYMinBoundingBox
The minimum y value of the bounding box.
Definition simulation.h:500
CSimulation::FindAllSeaCells
void FindAllSeaCells(void)
Finds and flags all sea areas which have at least one cell at a grid edge (i.e. does not flag 'inland...
Definition locate_coast.cpp:78
CSimulation::LSmoothCoastRunningMean
CGeomLine LSmoothCoastRunningMean(CGeomLine *) const
Does running-mean smoothing of a CGeomLine coastline vector (is in external CRS coordinates)
Definition smooth_line.cpp:209
CSimulation::nLocateSeaAndCoasts
int nLocateSeaAndCoasts(int &)
First find all connected sea areas, then locate the vector coastline(s), then put these onto the rast...
Definition locate_coast.cpp:55
CSimulation::m_dThisIterSWL
double m_dThisIterSWL
The still water level for this timestep (this includes tidal changes and any long-term SWL change)
Definition simulation.h:667
CSimulation::m_pRasterGrid
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
Definition simulation.h:1478
CSimulation::m_nXGridSize
int m_nXGridSize
The size of the grid in the x direction.
Definition simulation.h:431
CSimulation::nGetOppositeDirection
static int nGetOppositeDirection(int const)
Returns the opposite direction.
Definition gis_utils.cpp:1688
CSimulation::LogStream
ofstream LogStream
Definition simulation.h:1768
CSimulation::m_VFloodWaveSetupSurgeRunup
vector< CRWCoast > m_VFloodWaveSetupSurgeRunup
TODO 007 Info needed.
Definition simulation.h:1487
CSimulation::m_VCoast
vector< CRWCoast > m_VCoast
The coastline objects.
Definition simulation.h:1481
CSimulation::nLocateFloodAndCoasts
int nLocateFloodAndCoasts(void)
First find all connected sea areas, then locate the vector coastline(s), then put these onto the rast...
Definition locate_coast.cpp:1029
CSimulation::nTraceCoastLine
int nTraceCoastLine(unsigned int const, int const, int const, vector< bool > *, vector< CGeom2DIPoint > const *)
Traces a coastline (which is defined to be just above still water level) on the grid using the 'wall ...
Definition locate_coast.cpp:393
CSimulation::nTraceAllFloodCoasts
int nTraceAllFloodCoasts(void)
Locates all the potential coastline start points on the edges of the raster grid, then traces vector ...
Definition locate_coast.cpp:1265
CSimulation::FloodFillLand
void FloodFillLand(int const, int const)
Use the sealevel, wave set-up and run-up to evaluate flood hydraulically connected TODO 007 Not clear...
Definition locate_coast.cpp:1111
CSimulation::m_dThisIterDiffWaveSetupSurgeWaterLevel
double m_dThisIterDiffWaveSetupSurgeWaterLevel
TODO 007 Info needed.
Definition simulation.h:688
CSimulation::m_nYGridSize
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:434
CSimulation::m_dThisIterTopElevMin
double m_dThisIterTopElevMin
This-iteration lowest elevation of DEM.
Definition simulation.h:919
CSimulation::CalcCoastTangents
void CalcCoastTangents(int const)
Calculates tangents to a coastline: the tangent is assumed to be the orientation of energy/sediment f...
Definition calc_waves.cpp:2091
CSimulation::m_VFloodWaveSetupSurge
vector< CRWCoast > m_VFloodWaveSetupSurge
TODO 007 Info needed.
Definition simulation.h:1484
CSimulation::m_nCoastMax
int m_nCoastMax
Maximum valid coast length when searching for coasts, actually is COAST_LENGTH_MAX * tMax(m_nXGridSiz...
Definition simulation.h:470
CSimulation::m_nXMaxBoundingBox
int m_nXMaxBoundingBox
The maximum x value of the bounding box.
Definition simulation.h:497
CSimulation::m_VdFloodLocationX
vector< double > m_VdFloodLocationX
X coordinate (grid CRS) for total water level flooding.
Definition simulation.h:1382
CSimulation::m_bOmitSearchWestEdge
bool m_bOmitSearchWestEdge
Omit the west edge of the grid from coast-end searches?
Definition simulation.h:332
CSimulation::nTraceFloodCoastLine
int nTraceFloodCoastLine(unsigned int const, int const, int const, vector< bool > *, vector< CGeom2DIPoint > const *)
Traces a coastline (which is defined to be just above still water level) on the grid using the 'wall ...
Definition locate_coast.cpp:1362
CSimulation::dGridCentroidYToExtCRSY
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:71
CSimulation::m_nCoastMin
int m_nCoastMin
Minimum valid coast legth when searching for coass, actualli is tMin(m_nXGridSize,...
Definition simulation.h:473
CSimulation::m_bOmitSearchNorthEdge
bool m_bOmitSearchNorthEdge
Omit the north edge of the grid from coast-end searches?
Definition simulation.h:326
CSimulation::m_VEdgeCell
vector< CGeom2DIPoint > m_VEdgeCell
Edge cells.
Definition simulation.h:1493
CSimulation::m_nYMaxBoundingBox
int m_nYMaxBoundingBox
The maximum y value of the bounding box.
Definition simulation.h:503
CSimulation::m_ulThisIterNumSeaCells
unsigned long m_ulThisIterNumSeaCells
The number of grid cells which are marked as sea, for this iteration.
Definition simulation.h:565
CSimulation::LSmoothCoastSavitzkyGolay
CGeomLine LSmoothCoastSavitzkyGolay(CGeomLine *, int const, int const) const
Does smoothing of a CGeomLine coastline vector (is in external CRS coordinates) using a Savitzky-Gola...
Definition smooth_line.cpp:76
CSimulation::m_nXMinBoundingBox
int m_nXMinBoundingBox
The minimum x value of the bounding box.
Definition simulation.h:494
CSimulation::bIsWithinValidGrid
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,...
Definition gis_utils.cpp:170
CSimulation::m_nCoastSmooth
int m_nCoastSmooth
Which method to use for coast smoothing.
Definition simulation.h:440
CSimulation::DoCoastCurvature
void DoCoastCurvature(int const, int const)
Calculates both detailed and smoothed curvature for every point on a coastline.
Definition calc_curvature.cpp:43
CSimulation::m_bOmitSearchSouthEdge
bool m_bOmitSearchSouthEdge
Omit the south edge of the grid from coast-end searches?
Definition simulation.h:329
CSimulation::m_ulIter
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:553
CSimulation::dGridCentroidXToExtCRSX
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
CSimulation::nTraceAllCoasts
int nTraceAllCoasts(int &)
Locates all the potential coastline start points on the edges of the raster grid, then traces vector ...
Definition locate_coast.cpp:287
CSimulation::m_dThisIterDiffWaveSetupSurgeRunupWaterLevel
double m_dThisIterDiffWaveSetupSurgeRunupWaterLevel
TODO 007 Info needed.
Definition simulation.h:691
CSimulation::bIsInterventionCell
bool bIsInterventionCell(int const, int const) const
Returns true if the cell is an intervention.
Definition utils.cpp:2736
CSimulation::m_dThisIterTopElevMax
double m_dThisIterTopElevMax
This-iteration highest elevation of DEM.
Definition simulation.h:916
CSimulation::m_VdFloodLocationY
vector< double > m_VdFloodLocationY
X coordinate (grid CRS) for total water level flooding.
Definition simulation.h:1385
CSimulation::FloodFillSea
void FloodFillSea(int const, int const)
Flood-fills all sea cells starting from a given cell. The flood fill code used here is adapted from a...
Definition locate_coast.cpp:109
CSimulation::m_dThisIterDiffTotWaterLevel
double m_dThisIterDiffTotWaterLevel
TODO 007 Info needed.
Definition simulation.h:682
CSimulation::FindAllInundatedCells
int FindAllInundatedCells(void)
Finds and flags all sea areas which have at least one cell at a grid edge (i.e. does not flag 'inland...
Definition locate_coast.cpp:1069
CSimulation::m_nLevel
int m_nLevel
TODO 007 Used in WAVESETUP + SURGE + RUNUP.
Definition simulation.h:536
CSimulation::m_bOmitSearchEastEdge
bool m_bOmitSearchEastEdge
Omit the east edge of the grid from coast-end searches?
Definition simulation.h:335
CSimulation::m_VEdgeCellEdge
vector< int > m_VEdgeCellEdge
The grid edge that each edge cell belongs to.
Definition simulation.h:1496
cme.h
This file contains global definitions for CoastalME.
TOLERANCE
double const TOLERANCE
Definition cme.h:698
LF_CAT_SEA
int const LF_CAT_SEA
Definition cme.h:422
RTN_ERR_NOCOAST
int const RTN_ERR_NOCOAST
Definition cme.h:611
LF_CAT_SEDIMENT_INPUT
int const LF_CAT_SEDIMENT_INPUT
Definition cme.h:424
ERR
string const ERR
Definition cme.h:775
LF_CAT_SEDIMENT_INPUT_SUBMERGED
int const LF_CAT_SEDIMENT_INPUT_SUBMERGED
Definition cme.h:425
LOG_FILE_MIDDLE_DETAIL
int const LOG_FILE_MIDDLE_DETAIL
Definition cme.h:377
bFPIsEqual
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1176
LOG_FILE_HIGH_DETAIL
int const LOG_FILE_HIGH_DETAIL
Definition cme.h:378
SMOOTH_SAVITZKY_GOLAY
int const SMOOTH_SAVITZKY_GOLAY
Definition cme.h:657
SMOOTH_RUNNING_MEAN
int const SMOOTH_RUNNING_MEAN
Definition cme.h:656
SOUTH
int const SOUTH
Definition cme.h:387
LOG_FILE_ALL
int const LOG_FILE_ALL
Definition cme.h:379
EAST
int const EAST
Definition cme.h:385
RTN_OK
int const RTN_OK
Definition cme.h:577
LF_CAT_SEDIMENT_INPUT_NOT_SUBMERGED
int const LF_CAT_SEDIMENT_INPUT_NOT_SUBMERGED
Definition cme.h:426
RIGHT_HANDED
int const RIGHT_HANDED
Definition cme.h:397
RTN_ERR_TRACING_COAST
int const RTN_ERR_TRACING_COAST
Definition cme.h:610
tAbs
T tAbs(T a)
Definition cme.h:1148
NORTH
int const NORTH
Definition cme.h:383
LEFT_HANDED
int const LEFT_HANDED
Definition cme.h:398
WEST
int const WEST
Definition cme.h:389
coast.h
Contains CRWCoast definitions.
i_line.h
Contains CGeomILine definitions.
line.h
Contains CGeomLine definitions.
raster_grid.h
Contains CGeomRasterGrid definitions.
simulation.h
Contains CSimulation definitions.