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