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
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#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//===============================================================================================================================
56{
57 // Find all connected sea cells
59
60 // Find every coastline on the raster grid, mark raster cells, then create the vector coastline
61 int nRet = nTraceAllCoasts(nValidCoast);
62
63 if (nRet != RTN_OK)
64 return nRet;
65
66 // Have we created any coasts?
67 if (m_VCoast.empty())
68 {
69 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;
70 return RTN_ERR_NOCOAST;
71 }
72
73 return RTN_OK;
74}
75
76//===============================================================================================================================
78//===============================================================================================================================
80{
81 // Go along the list of edge cells
82 for (unsigned int n = 0; n < m_VEdgeCell.size(); n++)
83 {
85 continue;
86
88 continue;
89
91 continue;
92
94 continue;
95
96 int nX = m_VEdgeCell[n].nGetX();
97 int nY = m_VEdgeCell[n].nGetY();
98
99 if ((m_pRasterGrid->m_Cell[nX][nY].bIsInundated()) && (bFPIsEqual(m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth(), 0.0, TOLERANCE)))
100 // This edge cell is below SWL and sea depth remains set to zero
101 CellByCellFillSea(nX, nY);
102 }
103}
104
105//===============================================================================================================================
107//===============================================================================================================================
108void CSimulation::CellByCellFillSea(int const nXStart, int const nYStart)
109{
110 // For safety check
111 int nRoundLoopMax = m_nXGridSize * m_nYGridSize;
112
113 // Create an empty stack
114 stack<CGeom2DIPoint> PtiStack;
115
116 // Start at the given edge cell, push this onto the stack
117 PtiStack.push(CGeom2DIPoint(nXStart, nYStart));
118
119 // Then do the cell-by-cell fill loop until there are no more cell coordinates on the stack
120 int nRoundLoop = 0;
121
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?
152 {
153 if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
154 {
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
162 else
163 {
165 }
166 }
167
168 else
169 {
170 // No sediment input here, just mark as sea
171 m_pRasterGrid->m_Cell[nX][nY].SetInContiguousSea();
172 pLandform->SetLFCategory(LF_CAT_SEA);
173
174 // 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
175 m_pRasterGrid->m_Cell[nX][nY].SetWaveValuesToDeepWaterWaveValues();
176 }
177
178 // Now sort out the x-y extremities of the contiguous sea for the bounding box (used later in wave propagation)
179 if (nX < m_nXMinBoundingBox)
181
182 if (nX > m_nXMaxBoundingBox)
184
185 if (nY < m_nYMinBoundingBox)
187
188 if (nY > m_nYMaxBoundingBox)
190
191 // Update count
193
194 if ((! bSpanAbove) && (nY > 0) && (! m_pRasterGrid->m_Cell[nX][nY - 1].bBasementElevIsMissingValue()) && (m_pRasterGrid->m_Cell[nX][nY - 1].bIsInundated()))
195 {
196 PtiStack.push(CGeom2DIPoint(nX, nY - 1));
197 bSpanAbove = true;
198 }
199
200 else if (bSpanAbove && (nY > 0) && (! m_pRasterGrid->m_Cell[nX][nY - 1].bBasementElevIsMissingValue()) && (! m_pRasterGrid->m_Cell[nX][nY - 1].bIsInundated()))
201 {
202 bSpanAbove = false;
203 }
204
205 if ((! bSpanBelow) && (nY < m_nYGridSize - 1) && (! m_pRasterGrid->m_Cell[nX][nY + 1].bBasementElevIsMissingValue()) && (m_pRasterGrid->m_Cell[nX][nY + 1].bIsInundated()))
206 {
207 PtiStack.push(CGeom2DIPoint(nX, nY + 1));
208 bSpanBelow = true;
209 }
210
211 else if (bSpanBelow && (nY < m_nYGridSize - 1) && (! m_pRasterGrid->m_Cell[nX][nY + 1].bBasementElevIsMissingValue()) && (! m_pRasterGrid->m_Cell[nX][nY + 1].bIsInundated()))
212 {
213 bSpanBelow = false;
214 }
215
216 nX++;
217 }
218 }
219
220 // DEBUG CODE ===========================================================================================================
221 string strOutFile = m_strOutPath + "is_contiguous_sea_";
222 strOutFile += to_string(m_ulIter);
223 strOutFile += ".tif";
224
225 GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
226 GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
227 pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
228 pDataSet->SetGeoTransform(m_dGeoTransform);
229 double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
230 int n = 0;
231 for (int nY = 0; nY < m_nYGridSize; nY++)
232 {
233 for (int nX = 0; nX < m_nXGridSize; nX++)
234 {
235 pdRaster[n++] = m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea();
236 }
237 }
238
239 GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
240 pBand->SetNoDataValue(m_dMissingValue);
241 int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
242 if (nRet == CE_Failure)
243 return;
244
245 GDALClose(pDataSet);
246 delete[] pdRaster;
247 // DEBUG CODE ===========================================================================================================
248
249 // // DEBUG CODE ===========================================================================================================
250 // string strOutFile = m_strOutPath + "is_inundated_";
251 // strOutFile += to_string(m_ulIter);
252 // strOutFile += ".tif";
253 //
254 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
255 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
256 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
257 // pDataSet->SetGeoTransform(m_dGeoTransform);
258 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
259 //
260 // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
261 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
262 // pDataSet->SetGeoTransform(m_dGeoTransform);
263 //
264 // pdRaster = new double[m_nXGridSize * m_nYGridSize];
265 // int n = 0;
266 // for (int nY = 0; nY < m_nYGridSize; nY++)
267 // {
268 // for (int nX = 0; nX < m_nXGridSize; nX++)
269 // {
270 // pdRaster[n++] = m_pRasterGrid->m_Cell[nX][nY].bIsInundated();
271 // }
272 // }
273 //
274 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
275 // pBand = pDataSet->GetRasterBand(1);
276 // pBand->SetNoDataValue(m_dMissingValue);
277 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
278 // if (nRet == CE_Failure)
279 // return;
280 //
281 // GDALClose(pDataSet);
282 // delete[] pdRaster;
283 // // DEBUG CODE ===========================================================================================================
284
285 // // DEBUG CODE ===========================================================================================================
286 // LogStream << m_ulIter << ": cell-by-cell 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;
287
288 // LogStream << " m_nXMinBoundingBox = " << m_nXMinBoundingBox << " m_nXMaxBoundingBox = " << m_nXMaxBoundingBox << " m_nYMinBoundingBox = " << m_nYMinBoundingBox << " m_nYMaxBoundingBox = " << m_nYMaxBoundingBox << endl;
289 // // DEBUG CODE ===========================================================================================================
290}
291
292//===============================================================================================================================
294//===============================================================================================================================
295int CSimulation::nTraceAllCoasts(int &nValidCoast)
296{
298 LogStream << m_ulIter << ": Tracing coast" << endl;
299
300 vector<bool> VbPossibleStartCellLHEdge;
301 vector<bool> VbTraced;
302 vector<int> VnSearchDirection;
303 vector<CGeom2DIPoint> V2DIPossibleStartCell;
304
305 // Go along the list of edge cells and look for possible coastline start/finish cells
306 for (unsigned int n = 0; n < m_VEdgeCell.size() - 1; n++)
307 {
309 continue;
310
312 continue;
313
315 continue;
316
318 continue;
319
320 int nXThis = m_VEdgeCell[n].nGetX();
321 int nYThis = m_VEdgeCell[n].nGetY();
322 int nXNext = m_VEdgeCell[n + 1].nGetX();
323 int nYNext = m_VEdgeCell[n + 1].nGetY();
324
325 // Get "Is it sea?" information for 'this' and 'next' cells
326 bool bThisCellIsSea = m_pRasterGrid->m_Cell[nXThis][nYThis].bIsInContiguousSea();
327 bool bNextCellIsSea = m_pRasterGrid->m_Cell[nXNext][nYNext].bIsInContiguousSea();
328
329 // Are we at a coast?
330 if ((! bThisCellIsSea) && bNextCellIsSea)
331 {
332 // '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)?
333 if (! m_pRasterGrid->m_Cell[nXThis][nYThis].bIsPossibleCoastStartCell())
334 {
335 // It has not, so flag it
336 m_pRasterGrid->m_Cell[nXThis][nYThis].SetPossibleCoastStartCell();
337
339 LogStream << m_ulIter << ": Flagging [" << nXThis << "][" << nYThis << "] = {" << dGridCentroidXToExtCRSX(nXThis) << ", " << dGridCentroidYToExtCRSY(nYThis) << "} as possible coast start cell (left_handed edge)" << endl;
340
341 // And save it
342 V2DIPossibleStartCell.push_back(CGeom2DIPoint(nXThis, nYThis));
343 VbPossibleStartCellLHEdge.push_back(true);
344 VnSearchDirection.push_back(nGetOppositeDirection(m_VEdgeCellEdge[n]));
345 VbTraced.push_back(false);
346 }
347 }
348
349 else if (bThisCellIsSea && (! bNextCellIsSea))
350 {
351 // 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)?
352 if (! m_pRasterGrid->m_Cell[nXNext][nYNext].bIsPossibleCoastStartCell())
353 {
354 // It has not, so flag it
355 m_pRasterGrid->m_Cell[nXNext][nYNext].SetPossibleCoastStartCell();
356
358 LogStream << m_ulIter << ": Flagging [" << nXNext << "][" << nYNext << "] = {" << dGridCentroidXToExtCRSX(nXNext) << ", " << dGridCentroidYToExtCRSY(nYNext) << "} as possible coast start cell (right_handed edge)" << endl;
359
360 // And save it
361 V2DIPossibleStartCell.push_back(CGeom2DIPoint(nXNext, nYNext));
362 VbPossibleStartCellLHEdge.push_back(false);
363 VnSearchDirection.push_back(nGetOppositeDirection(m_VEdgeCellEdge[n + 1]));
364 VbTraced.push_back(false);
365 }
366 }
367 }
368
369 // Any possible coastline start/finish cells found?
370 if (V2DIPossibleStartCell.size() == 0)
371 {
372 LogStream << m_ulIter << ": no coastline start/finish points found after grid edges searched.";
373
375 {
376 LogStream << " Note that the following grid edges were not searched: " << (m_bOmitSearchNorthEdge ? "N " : "") << (m_bOmitSearchSouthEdge ? "S " : "") << (m_bOmitSearchWestEdge ? "W " : "") << (m_bOmitSearchEastEdge ? "E " : "");
377 }
378
379 LogStream << endl;
380
382 }
383
384 // Some possible coastline start/finish points were found
385 LogStream << m_ulIter << ": " << V2DIPossibleStartCell.size() << " possible coastline start/finish points found" << endl;
386
387 // OK, so trace from each of these possible start/finish points
388 for (unsigned int n = 0; n < V2DIPossibleStartCell.size(); n++)
389 {
390 if (! VbTraced[n])
391 {
392 int nRet = 0;
393
394 if (VbPossibleStartCellLHEdge[n])
395 {
396 nRet = nTraceCoastLine(n, VnSearchDirection[n], LEFT_HANDED, &VbTraced, &V2DIPossibleStartCell);
397 }
398
399 else
400 {
401 nRet = nTraceCoastLine(n, VnSearchDirection[n], RIGHT_HANDED, &VbTraced, &V2DIPossibleStartCell);
402 }
403
404 if (nRet == RTN_OK)
405 {
406 // We have a valid coastline starting from this possible start cell
407 VbTraced[n] = true;
408 nValidCoast++;
409 }
410 }
411 }
412
413 if (nValidCoast > 0)
414 return RTN_OK;
415
416 else
417 {
418 cerr << m_ulIter << ": no valid coasts found, see " << m_strLogFile << " for more information" << endl;
420 }
421}
422
423//===============================================================================================================================
425//===============================================================================================================================
426int CSimulation::nTraceCoastLine(unsigned int const nTraceFromStartCellIndex, int const nStartSearchDirection, int const nHandedness, vector<bool>*pVbTraced, vector<CGeom2DIPoint> const* pV2DIPossibleStartCell)
427{
428 bool bHitStartCell = false;
429 bool bAtCoast = false;
430 bool bHasLeftStartEdge = false;
431 bool bTooLong = false;
432 bool bOffEdge = false;
433 bool bRepeating = false;
434
435 int nStartX = pV2DIPossibleStartCell->at(nTraceFromStartCellIndex).nGetX();
436 int nStartY = pV2DIPossibleStartCell->at(nTraceFromStartCellIndex).nGetY();
437 int nX = nStartX;
438 int nY = nStartY;
439 int nSearchDirection = nStartSearchDirection;
440 int nRoundLoop = -1;
441 // nThisLen = 0;
442 // nLastLen = 0,
443 // nPreLastLen = 0;
444
445 // Temporary coastline as integer points (grid CRS)
446 CGeomILine ILTempGridCRS;
447
448 // Mark the start cell as coast and add it to the vector object
449 m_pRasterGrid->m_Cell[nStartX][nStartY].SetAsCoastline(true);
450 CGeom2DIPoint PtiStart(nStartX, nStartY);
451 ILTempGridCRS.Append(&PtiStart);
452
453 // 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
454 do
455 {
456 // // DEBUG CODE ==============================================================================================================
457 // LogStream << "Now at [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
458 // LogStream << "ILTempGridCRS is now:" << endl;
459 // for (int n = 0; n < ILTempGridCRS.nGetSize(); n++)
460 // LogStream << "[" << ILTempGridCRS[n].nGetX() << "][" << ILTempGridCRS[n].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[n].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[n].nGetY()) << "}" << endl;
461 // LogStream << "=================" << endl;
462 // // DEBUG CODE ==============================================================================================================
463
464 // Safety check
465 if (++nRoundLoop > m_nCoastMax)
466 {
467 bTooLong = true;
468
469 // LogStream << m_ulIter << ": abandoning coastline tracing from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "}, exceeded maximum search length (" << m_nCoastMax << ")" << endl;
470
471 // for (int n = 0; n < ILTempGridCRS.nGetSize(); n++)
472 // LogStream << "[" << ILTempGridCRS[n].nGetX() << "][" << ILTempGridCRS[n].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[n].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[n].nGetY()) << "}" << endl;
473 // LogStream << endl;
474
475 break;
476 }
477
478 // Another safety check
479 if ((nRoundLoop > 10) && (ILTempGridCRS.nGetSize() < 2))
480 {
481 // 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
482 bRepeating = true;
483
484 break;
485 }
486
487 // OK so far: so have we left the start edge?
488 if (! bHasLeftStartEdge)
489 {
490 // We have not yet left the start edge
491 if (((nStartSearchDirection == SOUTH) && (nY > nStartY)) || ((nStartSearchDirection == NORTH) && (nY < nStartY)) ||
492 ((nStartSearchDirection == EAST) && (nX > nStartX)) || ((nStartSearchDirection == WEST) && (nX < nStartX)))
493 bHasLeftStartEdge = true;
494
495 // Flag this cell to ensure that it is not chosen as a coastline start cell later
496 m_pRasterGrid->m_Cell[nX][nY].SetPossibleCoastStartCell();
497 // LogStream << "Flagging [" << nX << "][" << nY << "] as possible coast start cell NOT YET LEFT EDGE" << endl;
498 }
499
500 // 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
501 // if (bHasLeftStartEdge && bAtCoast)
502 {
503 for (unsigned int nn = 0; nn < pVbTraced->size(); nn++)
504 {
505 if ((nn != nTraceFromStartCellIndex) && (! pVbTraced->at(nn)))
506 {
507 // LogStream << "[" << pV2DIPossibleStartCell->at(nn).nGetX() << "][" << pV2DIPossibleStartCell->at(nn).nGetY() << "]" << endl;
508
509 if (bAtCoast && (nX == pV2DIPossibleStartCell->at(nn).nGetX()) && (nY == pV2DIPossibleStartCell->at(nn).nGetY()))
510 {
512 LogStream << m_ulIter << ": Possible coastline found, traced from [" << nStartX << "][" << nStartY << "] and hit another coast start cell at [" << nX << "][" << nY << "]" << endl;
513
514 pVbTraced->at(nn) = true;
515 bHitStartCell = true;
516 break;
517 }
518 }
519 }
520
521 // LogStream << endl;
522 }
523
524 if (bHitStartCell)
525 break;
526
527 // OK now sort out the next iteration of the search
528 int nXSeaward = 0;
529 int nYSeaward = 0;
530 int nSeawardNewDirection = 0;
531 int nXStraightOn = 0;
532 int nYStraightOn = 0;
533 int nXAntiSeaward = 0;
534 int nYAntiSeaward = 0;
535 int nAntiSeawardNewDirection = 0;
536 int nXGoBack = 0;
537 int nYGoBack = 0;
538 int nGoBackNewDirection = 0;
539
540 CGeom2DIPoint Pti(nX, nY);
541
542 // Set up the variables
543 switch (nHandedness)
544 {
545 case RIGHT_HANDED:
546
547 // 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
548 switch (nSearchDirection)
549 {
550 case NORTH:
551 // The sea is towards the RHS (E) of the coast, so first try to go right (to the E)
552 nXSeaward = nX + 1;
553 nYSeaward = nY;
554 nSeawardNewDirection = EAST;
555
556 // If can't do this, try to go straight on (to the N)
557 nXStraightOn = nX;
558 nYStraightOn = nY - 1;
559
560 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (W)
561 nXAntiSeaward = nX - 1;
562 nYAntiSeaward = nY;
563 nAntiSeawardNewDirection = WEST;
564
565 // As a last resort, go back (to the S)
566 nXGoBack = nX;
567 nYGoBack = nY + 1;
568 nGoBackNewDirection = SOUTH;
569
570 break;
571
572 case EAST:
573 // The sea is towards the RHS (S) of the coast, so first try to go right (to the S)
574 nXSeaward = nX;
575 nYSeaward = nY + 1;
576 nSeawardNewDirection = SOUTH;
577
578 // If can't do this, try to go straight on (to the E)
579 nXStraightOn = nX + 1;
580 nYStraightOn = nY;
581
582 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (N)
583 nXAntiSeaward = nX;
584 nYAntiSeaward = nY - 1;
585 nAntiSeawardNewDirection = NORTH;
586
587 // As a last resort, go back (to the W)
588 nXGoBack = nX - 1;
589 nYGoBack = nY;
590 nGoBackNewDirection = WEST;
591
592 break;
593
594 case SOUTH:
595 // The sea is towards the RHS (W) of the coast, so first try to go right (to the W)
596 nXSeaward = nX - 1;
597 nYSeaward = nY;
598 nSeawardNewDirection = WEST;
599
600 // If can't do this, try to go straight on (to the S)
601 nXStraightOn = nX;
602 nYStraightOn = nY + 1;
603
604 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (E)
605 nXAntiSeaward = nX + 1;
606 nYAntiSeaward = nY;
607 nAntiSeawardNewDirection = EAST;
608
609 // As a last resort, go back (to the N)
610 nXGoBack = nX;
611 nYGoBack = nY - 1;
612 nGoBackNewDirection = NORTH;
613
614 break;
615
616 case WEST:
617 // The sea is towards the RHS (N) of the coast, so first try to go right (to the N)
618 nXSeaward = nX;
619 nYSeaward = nY - 1;
620 nSeawardNewDirection = NORTH;
621
622 // If can't do this, try to go straight on (to the W)
623 nXStraightOn = nX - 1;
624 nYStraightOn = nY;
625
626 // If can't do either of these, try to go anti-seaward i.e. towards the LHS (S)
627 nXAntiSeaward = nX;
628 nYAntiSeaward = nY + 1;
629 nAntiSeawardNewDirection = SOUTH;
630
631 // As a last resort, go back (to the E)
632 nXGoBack = nX + 1;
633 nYGoBack = nY;
634 nGoBackNewDirection = EAST;
635
636 break;
637 }
638
639 break;
640
641 case LEFT_HANDED:
642
643 // 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
644 switch (nSearchDirection)
645 {
646 case NORTH:
647 // The sea is towards the LHS (W) of the coast, so first try to go left (to the W)
648 nXSeaward = nX - 1;
649 nYSeaward = nY;
650 nSeawardNewDirection = WEST;
651
652 // If can't do this, try to go straight on (to the N)
653 nXStraightOn = nX;
654 nYStraightOn = nY - 1;
655
656 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (E)
657 nXAntiSeaward = nX + 1;
658 nYAntiSeaward = nY;
659 nAntiSeawardNewDirection = EAST;
660
661 // As a last resort, go back (to the S)
662 nXGoBack = nX;
663 nYGoBack = nY + 1;
664 nGoBackNewDirection = SOUTH;
665
666 break;
667
668 case EAST:
669 // The sea is towards the LHS (N) of the coast, so first try to go left (to the N)
670 nXSeaward = nX;
671 nYSeaward = nY - 1;
672 nSeawardNewDirection = NORTH;
673
674 // If can't do this, try to go straight on (to the E)
675 nXStraightOn = nX + 1;
676 nYStraightOn = nY;
677
678 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (S)
679 nXAntiSeaward = nX;
680 nYAntiSeaward = nY + 1;
681 nAntiSeawardNewDirection = SOUTH;
682
683 // As a last resort, go back (to the W)
684 nXGoBack = nX - 1;
685 nYGoBack = nY;
686 nGoBackNewDirection = WEST;
687
688 break;
689
690 case SOUTH:
691 // The sea is towards the LHS (E) of the coast, so first try to go left (to the E)
692 nXSeaward = nX + 1;
693 nYSeaward = nY;
694 nSeawardNewDirection = EAST;
695
696 // If can't do this, try to go straight on (to the S)
697 nXStraightOn = nX;
698 nYStraightOn = nY + 1;
699
700 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (W)
701 nXAntiSeaward = nX - 1;
702 nYAntiSeaward = nY;
703 nAntiSeawardNewDirection = WEST;
704
705 // As a last resort, go back (to the N)
706 nXGoBack = nX;
707 nYGoBack = nY - 1;
708 nGoBackNewDirection = NORTH;
709
710 break;
711
712 case WEST:
713 // The sea is towards the LHS (S) of the coast, so first try to go left (to the S)
714 nXSeaward = nX;
715 nYSeaward = nY + 1;
716 nSeawardNewDirection = SOUTH;
717
718 // If can't do this, try to go straight on (to the W)
719 nXStraightOn = nX - 1;
720 nYStraightOn = nY;
721
722 // If can't do either of these, try to go anti-seaward i.e. towards the RHS (N)
723 nXAntiSeaward = nX;
724 nYAntiSeaward = nY - 1;
725 nAntiSeawardNewDirection = NORTH;
726
727 // As a last resort, go back (to the E)
728 nXGoBack = nX + 1;
729 nYGoBack = nY;
730 nGoBackNewDirection = EAST;
731
732 break;
733 }
734
735 break;
736 }
737
738 // 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?
739 if (bIsWithinValidGrid(nXSeaward, nYSeaward))
740 {
741 // It is, so check if the cell in the seaward direction is a sea cell
742 if (m_pRasterGrid->m_Cell[nXSeaward][nYSeaward].bIsInContiguousSea())
743 {
744 // There is sea in this seaward direction, so we are on the coast
745 bAtCoast = true;
746
747 // Has the current cell already marked been marked as a coast cell?
748 if (! m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
749 {
750 // Not already marked, is this an intervention cell with the top above SWL?
751 if ((bIsInterventionCell(nX, nY)) && (m_pRasterGrid->m_Cell[nX][nY].dGetInterventionTopElev() >= m_dThisIterSWL))
752 {
753 // It is, so mark as coast and add it to the vector object
754 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
755 ILTempGridCRS.Append(&Pti);
756 }
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
767 else
768 {
769 // The seaward cell is not a sea cell, so we will move to it next time
770 nX = nXSeaward;
771 nY = nYSeaward;
772
773 // And set a new search direction, to keep turning seaward
774 nSearchDirection = nSeawardNewDirection;
775 continue;
776 }
777 }
778
779 // 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?
780 if (bIsWithinValidGrid(nXStraightOn, nYStraightOn))
781 {
782 // It is, so check if there is sea immediately in front
783 if (m_pRasterGrid->m_Cell[nXStraightOn][nYStraightOn].bIsInContiguousSea())
784 {
785 // Sea is in front, so we are on the coast
786 bAtCoast = true;
787
788 // Has the current cell already marked been marked as a coast cell?
789 if (! m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
790 {
791 // Not already marked, is this an intervention cell with the top above SWL?
792 if ((bIsInterventionCell(nX, nY)) && (m_pRasterGrid->m_Cell[nX][nY].dGetInterventionTopElev() >= m_dThisIterSWL))
793 {
794 // It is, so mark as coast and add it to the vector object
795 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
796 ILTempGridCRS.Append(&Pti);
797 }
798
799 else if (m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() >= m_dThisIterSWL)
800 {
801 // The sediment top is above SWL so mark as coast and add it to the vector object
802 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
803 ILTempGridCRS.Append(&Pti);
804 }
805 }
806 }
807
808 else
809 {
810 // The straight-ahead cell is not a sea cell, so we will move to it next time
811 nX = nXStraightOn;
812 nY = nYStraightOn;
813
814 // The search direction remains unchanged
815 continue;
816 }
817 }
818
819 // 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?
820 if (bIsWithinValidGrid(nXAntiSeaward, nYAntiSeaward))
821 {
822 // It is, so check if there is sea in this anti-seaward cell
823 if (m_pRasterGrid->m_Cell[nXAntiSeaward][nYAntiSeaward].bIsInContiguousSea())
824 {
825 // There is sea on the anti-seaward side, so we are on the coast
826 bAtCoast = true;
827
828 // Has the current cell already marked been marked as a coast cell?
829 if (! m_pRasterGrid->m_Cell[nX][nY].bIsCoastline())
830 {
831 // Not already marked, is this an intervention cell with the top above SWL?
832 if ((bIsInterventionCell(nX, nY)) && (m_pRasterGrid->m_Cell[nX][nY].dGetInterventionTopElev() >= m_dThisIterSWL))
833 {
834 // It is, so mark as coast and add it to the vector object
835 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
836 ILTempGridCRS.Append(&Pti);
837 }
838
839 else if (m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() >= m_dThisIterSWL)
840 {
841 // The sediment top is above SWL so mark as coast and add it to the vector object
842 m_pRasterGrid->m_Cell[nX][nY].SetAsCoastline(true);
843 ILTempGridCRS.Append(&Pti);
844 }
845 }
846 }
847
848 else
849 {
850 // The anti-seaward cell is not a sea cell, so we will move to it next time
851 nX = nXAntiSeaward;
852 nY = nYAntiSeaward;
853
854 // And set a new search direction, to keep turning seaward
855 nSearchDirection = nAntiSeawardNewDirection;
856 continue;
857 }
858 }
859
860 // 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
861 if (bIsWithinValidGrid(nXGoBack, nYGoBack))
862 {
863 nX = nXGoBack;
864 nY = nYGoBack;
865
866 // And change the search direction
867 nSearchDirection = nGoBackNewDirection;
868 }
869
870 else
871 {
872 // Our final choice is not a valid cell, so give up
873 bOffEdge = true;
874 break;
875 }
876 } while (true);
877
878 // OK, we have finished tracing this coastline on the grid. But is the coastline too long or too short?
879 int nCoastSize = ILTempGridCRS.nGetSize();
880
881 if (bOffEdge)
882 {
884 LogStream << m_ulIter << ": Ignoring possible 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;
885
886 // Unmark these cells as coast cells
887 for (int n = 0; n < nCoastSize; n++)
888 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
889
891 }
892
893 if (bTooLong)
894 {
895 // Around loop too many times, so abandon this coastline
897 {
898 LogStream << m_ulIter << ": Ignoring possible coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since round loop " << nRoundLoop << " times (m_nCoastMax = " << m_nCoastMax << "), coastline size is " << nCoastSize;
899
900 if (nCoastSize > 0)
901 LogStream << ", it ended at [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "}";
902
903 LogStream << endl;
904 }
905
906 // Unmark these cells as coast cells
907 for (int n = 0; n < nCoastSize; n++)
908 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
909
911 }
912
913 if (bRepeating)
914 {
916 {
917 LogStream << m_ulIter << ": Ignoring possible coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "} since repeating, coastline size is " << nCoastSize;
918
919 if (nCoastSize > 0)
920 LogStream << ", it ended at [" << ILTempGridCRS[nCoastSize - 1].nGetX() << "][" << ILTempGridCRS[nCoastSize - 1].nGetY() << "] = {" << dGridCentroidXToExtCRSX(ILTempGridCRS[nCoastSize - 1].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[nCoastSize - 1].nGetY()) << "}";
921
922 LogStream << endl;
923 }
924
925 // Unmark these cells as coast cells
926 for (int n = 0; n < nCoastSize; n++)
927 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
928
930 }
931
932 if (nCoastSize == 0)
933 {
934 // Zero-length coastline, so abandon it
936 LogStream << m_ulIter << ": abandoning zero-length coastline from [" << nStartX << "][" << nStartY << "] = {" << dGridCentroidXToExtCRSX(nStartX) << ", " << dGridCentroidYToExtCRSY(nStartY) << "}" << endl;
937
939 }
940
941 if (nCoastSize < m_nCoastMin)
942 {
943 // The vector coastline is too small, so abandon it
945 LogStream << m_ulIter << ": Ignoring possible 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;
946
947 // Unmark these cells as coast cells
948 for (int n = 0; n < nCoastSize; n++)
949 m_pRasterGrid->m_Cell[ILTempGridCRS[n].nGetX()][ILTempGridCRS[n].nGetY()].SetAsCoastline(false);
950
952 }
953
954 // OK this new coastline is fine
955 int nEndX = nX;
956 int nEndY = nY;
957 int nCoastEndX = ILTempGridCRS[nCoastSize - 1].nGetX();
958 int nCoastEndY = ILTempGridCRS[nCoastSize - 1].nGetY();
959
960 if ((nCoastEndX != nEndX) || (nCoastEndY != nEndY))
961 {
962 // 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?
963 if (! m_pRasterGrid->m_Cell[nCoastEndX][nCoastEndY].bIsBoundingBoxEdge())
964 {
965 // The final cell in ILTempGridCRS is not a grid-edge cell, so add the grid-edge cell and mark the cell as coastline
966 ILTempGridCRS.Append(nEndX, nEndY);
967 nCoastSize++;
968
969 m_pRasterGrid->m_Cell[nEndX][nEndY].SetAsCoastline(true);
970 }
971 }
972
973 // Need to specify start edge and end edge for smoothing routines
974 int nStartEdge = m_pRasterGrid->m_Cell[nStartX][nStartY].nGetBoundingBoxEdge();
975 int nEndEdge = m_pRasterGrid->m_Cell[nEndX][nEndY].nGetBoundingBoxEdge();
976
977 // 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
978 CGeomLine LTempExtCRS;
979
980 for (int j = 0; j < nCoastSize; j++)
981 {
982 LTempExtCRS.Append(dGridCentroidXToExtCRSX(ILTempGridCRS[j].nGetX()), dGridCentroidYToExtCRSY(ILTempGridCRS[j].nGetY()));
983 }
984
985 // Now do some smoothing of the vector output, if desired
987 LTempExtCRS = LSmoothCoastRunningMean(&LTempExtCRS);
988
990 LTempExtCRS = LSmoothCoastSavitzkyGolay(&LTempExtCRS, nStartEdge, nEndEdge);
991
992 // // DEBUG CODE ==================================================================================================
993 // LogStream << "==================================" << endl;
994 // for (int j = 0; j < nCoastSize; j++)
995 // {
996 // LogStream << "{" << dGridCentroidXToExtCRSX(ILTempGridCRS[j].nGetX()) << ", " << dGridCentroidYToExtCRSY(ILTempGridCRS[j].nGetY()) << "}" << "\t{" << LTempExtCRS.dGetXAt(j) << ", " << LTempExtCRS.dGetYAt(j) << "}" << endl;
997 // }
998 // LogStream << "==================================" << endl;
999 // // DEBUG CODE ==================================================================================================
1000
1001 // Create a new coastline object and append to it the vector of coastline objects
1002 CRWCoast CoastTmp(this);
1003 m_VCoast.push_back(CoastTmp);
1004 int nCoast = static_cast<int>(m_VCoast.size()) - 1;
1005
1006 // Set the coastline (Ext CRS)
1007 m_VCoast[nCoast].SetCoastlineExtCRS(&LTempExtCRS);
1008
1009 // Set the coastline (Grid CRS)
1010 m_VCoast[nCoast].SetCoastlineGridCRS(&ILTempGridCRS);
1011
1012 // CGeom2DPoint PtLast(DBL_MIN, DBL_MIN);
1013 // for (int j = 0; j < nCoastSize; j++)
1014 // {
1015 // // Store the smoothed points (in external CRS) in the coast's m_LCoastlineExtCRS object, also append dummy values to the other attribute vectors
1016 // if (PtLast != &LTempExtCRS[j]) // Avoid duplicate points
1017 // {
1018 // m_VCoast[nCoast].AppendPointToCoastlineExtCRS(LTempExtCRS[j].dGetX(), LTempExtCRS[j].dGetY());
1019 //
1020 // // Also store the locations of the corresponding unsmoothed points (in raster grid CRS) in the coast's m_ILCellsMarkedAsCoastline vector
1021 // m_VCoast[nCoast].AppendCellMarkedAsCoastline(&ILTempGridCRS[j]);
1022 // }
1023 //
1024 // PtLast = LTempExtCRS[j];
1025 // }
1026
1027 // Set values for the coast's other attributes: set the coast's handedness, and start and end edges
1028 m_VCoast[nCoast].SetSeaHandedness(nHandedness);
1029 m_VCoast[nCoast].SetStartEdge(nStartEdge);
1030 m_VCoast[nCoast].SetEndEdge(nEndEdge);
1031
1033 {
1034 LogStream << m_ulIter << ": Valid 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;
1035
1036 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 ";
1037
1038 if (nStartEdge == NORTH)
1039 LogStream << "north";
1040
1041 else if (nStartEdge == SOUTH)
1042 LogStream << "south";
1043
1044 else if (nStartEdge == WEST)
1045 LogStream << "west";
1046
1047 else if (nStartEdge == EAST)
1048 LogStream << "east";
1049
1050 LogStream << " edge to the ";
1051
1052 if (nEndEdge == NORTH)
1053 LogStream << "north";
1054
1055 else if (nEndEdge == SOUTH)
1056 LogStream << "south";
1057
1058 else if (nEndEdge == WEST)
1059 LogStream << "west";
1060
1061 else if (nEndEdge == EAST)
1062 LogStream << "east";
1063
1064 LogStream << " edge" << endl;
1065 }
1066
1067 // LogStream << "-----------------" << endl;
1068 // for (int kk = 0; kk < m_VCoast.back().nGetCoastlineSize(); kk++)
1069 // 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;
1070 // LogStream << "-----------------" << endl;
1071
1072 // Next calculate the curvature of the vector coastline
1073 DoCoastCurvature(nCoast, nHandedness);
1074
1075 // Calculate values for the coast's flux orientation vector
1076 CalcCoastTangents(nCoast);
1077
1078 // And create the vector of pointers to coastline-normal objects
1079 m_VCoast[nCoast].CreateProfilesAtCoastPoints();
1080
1081 return RTN_OK;
1082}
1083
1084//===============================================================================================================================
1086//===============================================================================================================================
1088{
1089 // Find all connected sea cells
1091
1092 // Find every coastline on the raster grid, mark raster cells, then create the vector coastline
1093 int nRet = nTraceAllFloodCoasts();
1094
1095 if (nRet != RTN_OK)
1096 return nRet;
1097
1098 // Have we created any coasts?
1099 switch (m_nLevel)
1100 {
1101 case 0: // WAVESETUP + SURGE:
1102 {
1103 if (m_VFloodWaveSetupSurge.empty())
1104 {
1105 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;
1106 return RTN_ERR_NOCOAST;
1107 }
1108
1109 break;
1110 }
1111
1112 case 1: // WAVESETUP + SURGE + RUNUP:
1113 {
1114 if (m_VFloodWaveSetupSurgeRunup.empty())
1115 {
1116 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;
1117 return RTN_ERR_NOCOAST;
1118 }
1119
1120 break;
1121 }
1122 }
1123
1124 return RTN_OK;
1125}
1126
1127//===============================================================================================================================
1129//===============================================================================================================================
1131{
1132 for (int nX = 0; nX < m_nXGridSize; nX++)
1133 {
1134 for (int nY = 0; nY < m_nYGridSize; nY++)
1135 {
1136 m_pRasterGrid->m_Cell[nX][nY].UnSetCheckFloodCell();
1137 m_pRasterGrid->m_Cell[nX][nY].UnSetInContiguousFlood();
1138 m_pRasterGrid->m_Cell[nX][nY].SetAsFloodLine(false);
1139 }
1140 }
1141
1142 // Go along the list of edge cells
1143 for (unsigned int n = 0; n < m_VEdgeCell.size(); n++)
1144 {
1146 continue;
1147
1149 continue;
1150
1152 continue;
1153
1155 continue;
1156
1157 int nX = m_VEdgeCell[n].nGetX();
1158 int nY = m_VEdgeCell[n].nGetY();
1159
1160 if ((! m_pRasterGrid->m_Cell[nX][nY].bIsCellFloodCheck()) && (m_pRasterGrid->m_Cell[nX][nY].bIsInundated()))
1161 {
1162 // This edge cell is below SWL and sea depth remains set to zero
1163 FloodFillLand(nX, nY);
1164 }
1165 }
1166
1167 return RTN_OK;
1168}
1169
int nGetSize(void) const
Returns the number of integer point in the vector which represents this 2D shape.
Definition 2di_shape.cpp:70
void Append(CGeom2DIPoint const *)
Appends a new integer point to the vector which represents this 2D shape.
Definition 2di_shape.cpp:81
void Append(CGeom2DPoint const *)
Appends a point to this 2D shape.
Definition 2d_shape.cpp:67
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
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Definition 2di_point.cpp:42
Geometry class used to represent 2D vector integer line objects.
Definition i_line.h:31
Geometry class used to represent 2D vector line objects.
Definition line.h:31
Real-world class used to represent the landform of a cell.
int nGetLFCategory(void) const
Get the landform category.
void SetLFCategory(int const)
Set the landform category.
Real-world class used to represent coastline objects.
Definition coast.h:40
int m_nLogFileDetail
Definition simulation.h:591
int m_nYMinBoundingBox
The minimum y value of the bounding box.
Definition simulation.h:558
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...
CGeomLine LSmoothCoastRunningMean(CGeomLine *) const
Does running-mean smoothing of a CGeomLine coastline vector (is in external CRS coordinates)
int nLocateSeaAndCoasts(int &)
First find all connected sea areas, then locate the vector coastline(s), then put these onto the rast...
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
static int nGetOppositeDirection(int const)
Returns the opposite direction.
ofstream LogStream
vector< CRWCoast > m_VFloodWaveSetupSurgeRunup
TODO 007 Info needed.
vector< CRWCoast > m_VCoast
The coastline objects.
int nLocateFloodAndCoasts(void)
First find all connected sea areas, then locate the vector coastline(s), then put these onto the rast...
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 ...
int nTraceAllFloodCoasts(void)
Locates all the potential coastline start points on the edges of the raster grid, then traces vector ...
void FloodFillLand(int const, int const)
Use the sealevel, wave set-up and run-up to evaluate flood hydraulically connected TODO 007 Not clear...
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:465
double m_dThisIterTopElevMin
This-iteration lowest elevation of DEM.
double m_dGeoTransform[6]
Definition simulation.h:731
double m_dMissingValue
Used by CoastalME for floating-point missing values.
void CalcCoastTangents(int const)
Calculates tangents to a coastline: the tangent is assumed to be the orientation of energy/sediment f...
vector< CRWCoast > m_VFloodWaveSetupSurge
TODO 007 Info needed.
int m_nXMaxBoundingBox
The maximum x value of the bounding box.
Definition simulation.h:555
string m_strLogFile
Name of output log file.
bool m_bOmitSearchWestEdge
Omit the west edge of the grid from coast-end searches?
Definition simulation.h:355
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
bool m_bOmitSearchNorthEdge
Omit the north edge of the grid from coast-end searches?
Definition simulation.h:349
vector< CGeom2DIPoint > m_VEdgeCell
Edge cells.
string m_strOutPath
Path for all output files.
int m_nYMaxBoundingBox
The maximum y value of the bounding box.
Definition simulation.h:561
unsigned long m_ulThisIterNumSeaCells
The number of grid cells which are marked as sea, for this iteration.
Definition simulation.h:630
CGeomLine LSmoothCoastSavitzkyGolay(CGeomLine *, int const, int const) const
Does smoothing of a CGeomLine coastline vector (is in external CRS coordinates) using a Savitzky-Gola...
int m_nXMinBoundingBox
The minimum x value of the bounding box.
Definition simulation.h:552
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 m_nCoastSmooth
Which method to use for coast smoothing.
Definition simulation.h:471
void DoCoastCurvature(int const, int const)
Calculates both detailed and smoothed curvature for every point on a coastline.
bool m_bOmitSearchSouthEdge
Omit the south edge of the grid from coast-end searches?
Definition simulation.h:352
string m_strGDALBasementDEMProjection
GDAL projection string for the basement DEM raster file.
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
int nTraceAllCoasts(int &)
Locates all the potential coastline start/finish points on the edges of the raster grid,...
void CellByCellFillSea(int const, int const)
Cell-by-cell fills all sea cells starting from a given cell. The cell-by-cell fill (aka 'floodfill') ...
bool bIsInterventionCell(int const, int const) const
Returns true if the cell is an intervention.
Definition utils.cpp:2918
double m_dThisIterTopElevMax
This-iteration highest elevation of DEM.
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...
char ** m_papszGDALRasterOptions
Options for GDAL when handling raster files.
Definition simulation.h:456
int m_nLevel
TODO 007 Used in WAVESETUP + SURGE + RUNUP.
Definition simulation.h:597
bool m_bOmitSearchEastEdge
Omit the east edge of the grid from coast-end searches?
Definition simulation.h:358
vector< int > m_VEdgeCellEdge
The grid edge that each edge cell belongs to.
This file contains global definitions for CoastalME.
double const TOLERANCE
Definition cme.h:811
int const LF_CAT_SEA
Definition cme.h:534
int const RTN_ERR_NOCOAST
Definition cme.h:726
int const LF_CAT_SEDIMENT_INPUT
Definition cme.h:536
string const ERR
Definition cme.h:890
int const LF_CAT_SEDIMENT_INPUT_SUBMERGED
Definition cme.h:537
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
int const LOG_FILE_HIGH_DETAIL
Definition cme.h:490
int const SMOOTH_SAVITZKY_GOLAY
Definition cme.h:772
int const SMOOTH_RUNNING_MEAN
Definition cme.h:771
int const SOUTH
Definition cme.h:499
int const LOG_FILE_ALL
Definition cme.h:491
int const EAST
Definition cme.h:497
int const RTN_OK
Definition cme.h:692
int const LF_CAT_SEDIMENT_INPUT_NOT_SUBMERGED
Definition cme.h:538
int const RIGHT_HANDED
Definition cme.h:509
int const RTN_ERR_TRACING_COAST
Definition cme.h:725
int const NORTH
Definition cme.h:495
int const LEFT_HANDED
Definition cme.h:510
int const WEST
Definition cme.h:501
Contains CRWCoast definitions.
Contains CGeomILine definitions.
Contains CGeomLine definitions.
Contains CGeomRasterGrid definitions.
Contains CSimulation definitions.