CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
create_polygons.cpp
Go to the documentation of this file.
1
12
13/* ==============================================================================================================================
14
15 This file is part of CoastalME, the Coastal Modelling Environment.
16
17 CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
20
21 You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22
23==============================================================================================================================*/
24#include <assert.h>
25
26#include <iostream>
27using std::endl;
28
29#include <stack>
30using std::stack;
31
32#include "cme.h"
33#include "simulation.h"
34#include "coast.h"
35#include "2d_point.h"
36#include "2di_point.h"
37
38//===============================================================================================================================
40//===============================================================================================================================
42{
43 // Global polygon count TODO 044
45
46 // Do this for each coast
47 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
48 {
49 int nNode = -1;
50 int nNextProfile = -1;
51 int nPolygon = -1;
52
53 for (int i = 0; i < static_cast<int>(m_pVCoastPolygon.size()); i++)
54 delete m_pVCoastPolygon[i];
55
56 m_pVCoastPolygon.clear();
57
58 // Do this for every point on the coastline (except the last point)
59 int const nCoastSize = m_VCoast[nCoast].nGetCoastlineSize();
60
61 for (int nCoastPoint = 0; nCoastPoint < nCoastSize - 1; nCoastPoint++)
62 {
63 if (! m_VCoast[nCoast].bIsProfileAtCoastPoint(nCoastPoint))
64 continue;
65
66 // OK, this coast point is the start of a coastline-normal profile
67 CGeomProfile* pThisProfile = m_VCoast[nCoast].pGetProfileAtCoastPoint(nCoastPoint);
68
69 if (pThisProfile->bOKIncStartAndEndOfCoast())
70 {
71 // This profile is OK, so we will start a polygon here and extend it down-coast (i.e. along the coast in the direction of increasing coastline point numbers)
72 int const nThisProfile = pThisProfile->nGetProfileID();
73
74 // This will be the coast ID number of the polygon, and also the polygon's along-coast sequence
75 nPolygon++;
76
77 // Now get a pointer to the next (down-coast) profile
78 CGeomProfile* pNextProfile = pThisProfile->pGetDownCoastAdjacentProfile();
79
80 bool bNextProfileIsOK = false;
81
82 do
83 {
84 // if (pNextProfile == NULL) // TODO THIS GIVES CPPCHECK ERROR
85 // {
86 // LogStream << m_ulIter << ": nThisProfile = " << nThisProfile << " invalid down-coast adjacent profile, trying next down-coast adjacent profile" << endl;
87 //
88 // // Try the next-after-next profile
89 // CGeomProfile* pNextNextProfile = pNextProfile->pGetDownCoastAdjacentProfile();
90 // pNextProfile = pNextNextProfile;
91 // continue;
92 // }
93
94 // Get the ID of the next (down-coast) profile
95 nNextProfile = pNextProfile->nGetProfileID();
96
97 // Is the next profile OK?
98 bNextProfileIsOK = pNextProfile->bOKIncStartAndEndOfCoast();
99
100 if (! bNextProfileIsOK)
101 {
102 // Nope, the next profile is not OK
103 // LogStream << m_ulIter << ": down-coast adjacent profile = " << nNextProfile << " is not OK" << endl;
104
105 // So try the following profile
106 CGeomProfile* pNextNextProfile = pNextProfile->pGetDownCoastAdjacentProfile();
107 pNextProfile = pNextNextProfile;
108 }
109
110 } while (! bNextProfileIsOK);
111
112 // LogStream << "Profile " << pNextProfile->nGetProfileID() << " is OK" << endl;
113
114 // Get the coast point at which this next profile starts
115 int const nNextProfileCoastPoint = pNextProfile->nGetCoastPoint();
116
117 // Calculate half the along-coast distance (in coast points) between this profile and the next (i.e. down-coast) profile
118 int const nDist = (nNextProfileCoastPoint - nCoastPoint) / 2;
119
120 // OK, set the node point in the coast object. We do this now, instead of earlier on, since some profiles (i.e. polygon boundaries) may have been marked as invalid
121 int const nNodePoint = nCoastPoint + nDist;
122 m_VCoast[nCoast].SetPolygonNode(nNodePoint, ++nNode);
123
124 // Get the grid CRS coordinates of the coast node
125 CGeom2DIPoint const PtiNode = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nNodePoint);
126
127 // // DEBUG CODE ===================
128 // if ((m_ulIter == 18) && (nThisProfile == 29))
129 // std::cerr << endl;
130 // // DEBUG CODE ===================
131
132 // Get some defaults (assuming for now that this polygon is not approximately triangular i.e. both normals do not meet)
133 int nNextProfileEnd = pNextProfile->nGetProfileSize() - 1;
134 int nThisProfileEnd = pThisProfile->nGetProfileSize() - 1;
135 bool bMeetsAtAPoint = false;
136 CGeom2DPoint PtCoastwardTip;
137
138 // // DEBUG CODE =============================================================================================
139 // CGeom2DPoint PtNextEndTmp = *pNextProfile->pPtGetPointInProfile(nNextProfileEnd);
140 // double dXTmp = dExtCRSXToGridX(PtNextEndTmp.dGetX());
141 // double dYTmp = dExtCRSYToGridY(PtNextEndTmp.dGetY());
142 // assert(bIsWithinValidGrid(dXTmp, dYTmp));
143 //
144 // CGeom2DPoint PtThisEndTmp = *pThisProfile->pPtGetPointInProfile(nThisProfileEnd);
145 // dXTmp = dExtCRSXToGridX(PtThisEndTmp.dGetX());
146 // dYTmp = dExtCRSYToGridY(PtThisEndTmp.dGetY());
147 // assert(bIsWithinValidGrid(dXTmp, dYTmp));
148 // // DEBUG CODE =============================================================================================
149
150 // Now check to see if the two normals do meet i.e. if they are coincident
151 if (pThisProfile->bFindProfileInCoincidentProfiles(nNextProfile))
152 {
153 // Yes they do meet
154 bMeetsAtAPoint = true;
155 int nTmpThisProfileEnd;
156 int nTmpNextProfileEnd;
157
158 // Find the most coastward point at which this normal and the previous normal touch. If they do not touch, the polygon requires a 'joining line'
159 pThisProfile->GetMostCoastwardSharedLineSegment(nNextProfile, nTmpThisProfileEnd, nTmpNextProfileEnd);
160
161 if (nTmpThisProfileEnd == -1)
162 {
163 LogStream << m_ulIter << ": " << ERR << "profile " << nNextProfile << " should be coincident with profile " << nThisProfile << " but was not found" << endl;
165 }
166
167 // Safety check: make sure that nThisProfileEnd is no bigger than pThisProfile->nGetProfileSize()-1, and the same for nNextProfileEnd
168 nThisProfileEnd = tMin(nThisProfileEnd, nTmpThisProfileEnd);
169 nNextProfileEnd = tMin(nNextProfileEnd, nTmpNextProfileEnd);
170
171 PtCoastwardTip = *pThisProfile->pPtGetPointInProfile(nThisProfileEnd);
172 }
173
174 // Create the vector in which to store the polygon's boundary (external CRS). Note that these points are not in sequence
175 vector<CGeom2DPoint> PtVBoundary;
176
177 // Start appending points: begin by appending the points in this normal, in reverse order
178 for (int i = nThisProfileEnd; i >= 0; i--)
179 {
180 CGeom2DPoint const PtThis = *pThisProfile->pPtGetPointInProfile(i);
181 PtVBoundary.push_back(PtThis);
182 }
183
184 // Next add coast points: from the start point of this normal, moving down-coast as far as the down-coast norma
185 for (int i = nCoastPoint; i <= nNextProfileCoastPoint; i++)
186 {
187 CGeom2DPoint const PtThis = *m_VCoast[nCoast].pPtGetCoastlinePointExtCRS(i);
188 PtVBoundary.push_back(PtThis);
189 }
190
191 // Append the points in the down-coast (next) normal
192 for (int i = 0; i <= nNextProfileEnd; i++)
193 {
194 CGeom2DPoint const PtThis = *pNextProfile->pPtGetPointInProfile(i);
195 PtVBoundary.push_back(PtThis);
196 }
197
198 // Finally, append the end point of this normal
199 CGeom2DPoint const PtThis = *pThisProfile->pPtGetPointInProfile(nThisProfileEnd);
200 PtVBoundary.push_back(PtThis);
201
202 // Now identify the 'anti-node', this is the seaward point 'opposite' the polygon's coastal node
203 CGeom2DIPoint PtiAntiNode;
204
205 if (bMeetsAtAPoint)
206 PtiAntiNode = PtiExtCRSToGridRound(&PtCoastwardTip);
207
208 else
209 {
210 CGeom2DPoint const PtAvg = PtAverage(pThisProfile->pPtGetPointInProfile(nThisProfileEnd), pNextProfile->pPtGetPointInProfile(nNextProfileEnd));
211 PtiAntiNode = PtiExtCRSToGridRound(&PtAvg);
212 }
213
214 // Is this a start-of-coast or an end-of-coast polygon?
215 bool bStartCoast = false;
216 bool bEndCoast = false;
217
218 if (pThisProfile->bStartOfCoast())
219 bStartCoast = true;
220
221 if (pNextProfile->bEndOfCoast())
222 bEndCoast = true;
223
224 // Create the coast polygon object and get a pointer to it. TODO 044 the first parameter (global ID) will need to change when considering multiple coasts
225 CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pPolyCreatePolygon(nPolygon, nNodePoint, &PtiNode, &PtiAntiNode, nThisProfile, nNextProfile, &PtVBoundary, nThisProfileEnd + 1, nNextProfileEnd + 1, bStartCoast, bEndCoast);
226
227 // And store this pointer for simulation-wide access, in along-coast sequence
228 m_pVCoastPolygon.push_back(pPolygon);
229
230 // Save the profile-end vertices (but omit the last one if the profiles meet at a point)
231 pPolygon->AppendVertex(pThisProfile->pPtiGetStartPoint());
232 pPolygon->AppendVertex(pThisProfile->pPtiGetEndPoint());
233 pPolygon->AppendVertex(pNextProfile->pPtiGetStartPoint());
234
235 if (! bMeetsAtAPoint)
236 pPolygon->AppendVertex(pNextProfile->pPtiGetEndPoint());
237
238 // // DEBUG CODE =================================================================================
239 // LogStream << m_ulIter << ": vertices for polygon = " << nPolygon << " (m_nNumPolygonGlobal = " << m_nNumPolygonGlobal << ")" << endl;
240 // for (int n = 0; n < pPolygon->nGetNumVertices(); n++)
241 // LogStream << "[" << pPolygon->PtiGetVertex(n).nGetX() << "][" << pPolygon->PtiGetVertex(n).nGetY() << "]\t";
242 // LogStream << endl;
243 // // DEBUG CODE =================================================================================
244
245 // assert(nPolygon < m_VCoast[nCoast].nGetNumPolygons());
246
247 // Now rasterize the polygon boundaries: first, the coastline. This is necessary so that sand/coarse sediment derived from platform erosion of the coast cells is correctly added to the containing polygon's unconsolidated sediment
248 for (int i = nCoastPoint; i <= nNextProfileCoastPoint; i++)
249 {
250 CGeom2DIPoint const PtiToMark = *m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(i);
251 m_pRasterGrid->m_Cell[PtiToMark.nGetX()][PtiToMark.nGetY()].SetCoastAndPolygonID(nCoast, nPolygon);
252 }
253
254 // Do the down-coast normal profile (does the whole length, including any shared line segments. So some cells are marked twice, however this is not a problem)
255 int nCellsInProfile = pNextProfile->nGetNumCellsInProfile();
256
257 for (int i = 0; i < nCellsInProfile; i++)
258 {
259 CGeom2DIPoint const PtiToMark = *pNextProfile->pPtiGetCellInProfile(i);
260 m_pRasterGrid->m_Cell[PtiToMark.nGetX()][PtiToMark.nGetY()].SetCoastAndPolygonID(nCoast, nPolygon);
261 }
262
263 // Do this normal profile (again does the whole length, including any shared line segments)
264 nCellsInProfile = pThisProfile->nGetNumCellsInProfile();
265
266 for (int i = 0; i < nCellsInProfile; i++)
267 {
268 CGeom2DIPoint const PtiToMark = *pThisProfile->pPtiGetCellInProfile(i);
269 m_pRasterGrid->m_Cell[PtiToMark.nGetX()][PtiToMark.nGetY()].SetCoastAndPolygonID(nCoast, nPolygon);
270 }
271
272 // If the polygon doesn't meet at a point at its seaward end, also need to rasterize the 'joining line'
273 if (! bMeetsAtAPoint)
274 {
275 CGeom2DIPoint const PtiDownCoastNormalEnd = *pNextProfile->pPtiGetEndPoint(); // Grid CRS
276 CGeom2DIPoint const PtiUpCoastNormalEnd = *pThisProfile->pPtiGetEndPoint(); // Grid CRS
277
278 RasterizePolygonJoiningLine(nCoast, &PtiUpCoastNormalEnd, &PtiDownCoastNormalEnd, nPolygon);
279
280 // pPolygon->SetNotPointed();
281 }
282
283 // }
284
285 // nNextProfile = nThisProfile;
286 nCoastPoint = nNextProfileCoastPoint - 1;
287
289
290 if (bEndCoast)
291 break;
292 }
293 }
294 }
295
296 return RTN_OK;
297}
298
299//===============================================================================================================================
301//===============================================================================================================================
302void CSimulation::RasterizePolygonJoiningLine(int nCoast, CGeom2DIPoint const* pPt1, CGeom2DIPoint const* pPt2, int const nPoly)
303{
304 // The start point of the line (grid CRS)
305 int const nXStart = pPt1->nGetX();
306 int const nYStart = pPt1->nGetY();
307
308 // The end point of the line (grid CRS)
309 int const nXEnd = pPt2->nGetX();
310 int const nYEnd = pPt2->nGetY();
311
312 // Safety check, in case the two points are identical (can happen due to rounding errors)
313 if ((nXStart == nXEnd) && (nYStart == nYEnd))
314 return;
315
316 // Interpolate between cells by a simple DDA line algorithm, see http://en.wikipedia.org/wiki/Digital_differential_analyzer_(graphics_algorithm) Note that Bresenham's algorithm gave occasional gaps
317 double dXInc = nXEnd - nXStart;
318 double dYInc = nYEnd - nYStart;
319 double const dLength = tMax(tAbs(dXInc), tAbs(dYInc));
320
321 dXInc /= dLength;
322 dYInc /= dLength;
323
324 double dX = nXStart;
325 double dY = nYStart;
326
327 // Process each interpolated point
328 for (int m = 0; m <= nRound(dLength); m++)
329 {
330 int nX = nRound(dX);
331 int nY = nRound(dY);
332
333 // Safety check
334 if (! bIsWithinValidGrid(nX, nY))
335 KeepWithinValidGrid(nXStart, nYStart, nX, nY);
336
337 // assert(nPoly < m_VCoast[0].nGetNumPolygons());
338
339 // Mark this point on the raster grid
340 m_pRasterGrid->m_Cell[nX][nY].SetCoastAndPolygonID(nCoast, nPoly);
341
342 // And increment for next time
343 dX += dXInc;
344 dY += dYInc;
345 }
346}
347
348//===============================================================================================================================
350//===============================================================================================================================
352{
353 // // DEBUG CODE =================================================================================
354 // vector<CGeom2DIPoint> VPtiCentroids;
355 // vector<int> VnID;
356 //
357 // // Do this for each coast
358 // for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
359 // {
360 // // Do this for every coastal polygon
361 // for (int nPoly = 0; nPoly < m_VCoast[nCoast].nGetNumPolygons(); nPoly++)
362 // {
363 // CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
364 //
365 // int nPolyID = pPolygon->nGetPolygonCoastID();
366 // VnID.push_back(nPolyID);
367 //
368 // CGeom2DIPoint PtiStart = pPolygon->PtiGetFillStartPoint();
369 // VPtiCentroids.push_back(PtiStart);
370 // }
371 // }
372 //
373 // string strOutFile = m_strOutPath;
374 // strOutFile += "00_polygon_flood_fill_start_point_";
375 // strOutFile += to_string(m_ulIter);
376 // strOutFile += ".tif";
377 //
378 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
379 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
380 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
381 // pDataSet->SetGeoTransform(m_dGeoTransform);
382 //
383 // int nn = 0;
384 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
385 // for (int nY = 0; nY < m_nYGridSize; nY++)
386 // {
387 // for (int nX = 0; nX < m_nXGridSize; nX++)
388 // {
389 // pdRaster[nn] = INT_NODATA;
390 //
391 // for (int n = 0; n < static_cast<int>(VPtiCentroids.size()); n++)
392 // {
393 // if ((VPtiCentroids[n].nGetX() == nX) && (VPtiCentroids[n].nGetY() == nY))
394 // {
395 // pdRaster[nn] = VnID[n];
396 // LogStream << m_ulIter << ": cell-by-cell fill start point for polygon " << VnID[n] << " is [" << nX << "][" << nY << "]" << endl;
397 // break;
398 // }
399 // }
400 // nn++;
401 // }
402 // }
403 //
404 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
405 // pBand->SetNoDataValue(m_dMissingValue);
406 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
407 //
408 // if (nRet == CE_Failure)
409 // LogStream << nRet << endl;
410 //
411 // GDALClose(pDataSet);
412 // delete[] pdRaster;
413 // // DEBUG CODE ===========================================================================================================
414
415 // Do this for each coast
416 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
417 {
418 // Do this for every coastal polygon, in along-coast sequence
419 int const nPolygons = m_VCoast[nCoast].nGetNumPolygons();
420
421 for (int nPoly = 0; nPoly < nPolygons; nPoly++)
422 {
423 int nCellsInPolygon = 0;
424 double dTotDepth = 0;
425 double dStoredUnconsFine = 0;
426 double dStoredUnconsSand = 0;
427 double dStoredUnconsCoarse = 0;
428 double dStoredConsFine = 0;
429 double dStoredConsSand = 0;
430 double dStoredConsCoarse = 0;
431 double dSedimentInputFine = 0;
432 double dSedimentInputSand = 0;
433 double dSedimentInputCoarse = 0;
434
435 CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
436 int const nPolyID = pPolygon->nGetPolygonCoastID(); // TODO 044
437
438 // LogStream << m_ulIter << ": in MarkPolygonCells() nPoly = " << nPoly << " nPolyID = " << nPolyID << endl;
439
440 // Create an empty stack
441 stack<CGeom2DIPoint> PtiStack;
442
443 // // Since the polygon's vector boundary does not coincide exactly with the polygon's raster boundary, and the point-in-polygon check gives an indeterminate result if the point is exactly on the polygon's boundary, for safety we must construct a vector 'inner buffer' which is smaller than, and inside, the vector boundary TODO *** STILL NEEDED?
444 // int nHand = m_VCoast[nCoast].nGetSeaHandedness();
445 // int nSize = pPolygon->nGetBoundarySize();
446 // vector<CGeom2DPoint> PtVInnerBuffer;
447 //
448 // for (int i = 0; i < nSize-1; i++)
449 // {
450 // int j = i + 1;
451 // if (i == nSize-2) // We must ignore the duplicated node point
452 // j = 0;
453 //
454 // CGeom2DPoint PtThis = *pPolygon->pPtGetBoundaryPoint(i);
455 // CGeom2DPoint PtNext = *pPolygon->pPtGetBoundaryPoint(j);
456 //
457 // // Safety check
458 // if (PtThis == PtNext)
459 // continue;
460 //
461 // CGeom2DPoint PtBuffer = PtGetPerpendicular(&PtThis, &PtNext, m_dCellSide, nHand);
462 //
463 // PtVInnerBuffer.push_back(PtBuffer);
464 // }
465
466 // // DEBUG CODE ==============================================================================================
467 // LogStream << endl << m_ulIter << ": coast " << nCoast << ", polygon " << nPoly << endl;
468 // LogStream << "Boundary\t\t\tBuffer" << endl;
469 // for (int i = 0; i < pPolygon->nGetBoundarySize()-1; i++)
470 // LogStream << "{" << pPolygon->pPtGetBoundaryPoint(i)->dGetX() << ", " << pPolygon->pPtGetBoundaryPoint(i)->dGetY() << "}\t{" << PtVInnerBuffer[i].dGetX() << ", " << PtVInnerBuffer[i].dGetY() << "}" << endl;
471 // LogStream << endl;
472 // // DEBUG CODE ==============================================================================================
473
474 // Use the centroid as the start point for the cell-by-cell fill procedure
475 CGeom2DIPoint const PtiStart = pPolygon->PtiGetFillStartPoint();
476
477 // // Is the centroid within the inner buffer?
478 // if (! bIsWithinPolygon(&PtStart, &PtVInnerBuffer))
479 // {
480 // // No, it is not: the polygon must be a concave polygon. So keep looking for a point which is definitely inside the polygon, using an alternative method
481 // PtStart = PtFindPointInPolygon(&PtVInnerBuffer, pPolygon->nGetPointInPolygonSearchStartPoint());
482 // }
483 //
484 // // Safety check (PtFindPointInPolygon() returns CGeom2DPoint(DBL_NODATA, DBL_NODATA) if it cannot find a valid start point)
485 // if (bFPIsEqual(PtStart.dGetX(), DBL_NODATA, TOLERANCE))
486 // {
487 // LogStream << m_ulIter << ": " << ERR << "could not find a cell-by-cell fill start point for coast " << nCoast << ", polygon " << nPoly << endl;
488 // break;
489 // }
490
491 // We have a cell-by-cell fill start point which is definitely within the polygon so push this point onto the stack
492 // CGeom2DIPoint PtiStart;
493 // PtiStart.SetX(nRound(PtStart.dGetX())); // Grid CRS
494 // PtiStart.SetY(nRound(PtStart.dGetY())); // Grid CRS
495
496 PtiStack.push(PtiStart);
497
498 // LogStream << m_ulIter << ": filling polygon " << nPoly << " from [" << PtiStart.nGetX() << "][" << PtiStart.nGetY() << "] = {" << dGridCentroidXToExtCRSX(PtiStart.nGetX()) << ", " << dGridCentroidYToExtCRSY(PtiStart.nGetY()) << "}" << endl;
499
500 // Then do the cell-by-cell fill: loop until there are no more cell coordinates on the stack
501 while (! PtiStack.empty())
502 {
503 CGeom2DIPoint const Pti = PtiStack.top();
504 PtiStack.pop();
505
506 int nX = Pti.nGetX();
507 int const nY = Pti.nGetY();
508
509 while ((nX >= 0) && (m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID() == INT_NODATA))
510 nX--;
511
512 nX++;
513
514 bool bSpanAbove = false;
515 bool bSpanBelow = false;
516
517 while ((nX < m_nXGridSize) && (m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID() == INT_NODATA))
518 {
519 // assert(nPolyID < m_VCoast[nCoast].nGetNumPolygons());
520
521 // Mark the cell as being in this polygon and this coast
522 m_pRasterGrid->m_Cell[nX][nY].SetCoastAndPolygonID(nCoast, nPolyID);
523
524 // LogStream << "Marked [" << nX << "][" << nY << "] = {" << dGridCentroidXToExtCRSX(nX) << ", " << dGridCentroidYToExtCRSY(nY) << "}" << endl;
525
526 // Increment the running totals for this polygon
527 // dCliffCollapseErosionFine += m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseErosionFine();
528 // dCliffCollapseErosionSand += m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseErosionSand();
529 // dCliffCollapseErosionCoarse += m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseErosionCoarse();
530 // dCliffCollapseTalusSand += m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseSandTalusDeposition();
531 // dCliffCollapseTalusCoarse += m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseCoarseTalusDeposition();
532
533 // Get the number of the highest layer with non-zero thickness
534 int const nThisLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
535
536 // Safety check
537 if ((nThisLayer != NO_NONZERO_THICKNESS_LAYERS) && (nThisLayer != INT_NODATA))
538 {
539 // Increment some more running totals for this polygon TODO 066 should this be for ALL layers above the basement?
540 dStoredUnconsFine += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetFineDepth();
541 dStoredUnconsSand += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
542 dStoredUnconsCoarse += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
543
544 dStoredConsFine += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetFineDepth();
545 dStoredConsSand += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetSandDepth();
546 dStoredConsCoarse += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetConsolidatedSediment()->dGetCoarseDepth();
547
548 // Add to the start-iteration total of suspended fine sediment within polygons
549 m_dStartIterSuspFineInPolygons += m_pRasterGrid->m_Cell[nX][nY].dGetSuspendedSediment();
550
551 // Add to the total of sediment derived from sediment input events
552 dSedimentInputFine += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetFineSedimentInputDepth();
553 dSedimentInputSand += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetSandSedimentInputDepth();
554 dSedimentInputCoarse += m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nThisLayer)->pGetUnconsolidatedSediment()->dGetCoarseSedimentInputDepth();
555 }
556
557 nCellsInPolygon++;
558 dTotDepth += m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth();
559
560 if ((! bSpanAbove) && (nY > 0) && (m_pRasterGrid->m_Cell[nX][nY - 1].nGetPolygonID() == INT_NODATA))
561 {
562 PtiStack.push(CGeom2DIPoint(nX, nY - 1));
563 bSpanAbove = true;
564 }
565
566 else if (bSpanAbove && (nY > 0) && (m_pRasterGrid->m_Cell[nX][nY - 1].nGetPolygonID() != INT_NODATA))
567 {
568 bSpanAbove = false;
569 }
570
571 if ((! bSpanBelow) && (nY < m_nYGridSize - 1) && (m_pRasterGrid->m_Cell[nX][nY + 1].nGetPolygonID() == INT_NODATA))
572 {
573 PtiStack.push(CGeom2DIPoint(nX, nY + 1));
574 bSpanBelow = true;
575 }
576
577 else if (bSpanBelow && (nY < m_nYGridSize - 1) && (m_pRasterGrid->m_Cell[nX][nY + 1].nGetPolygonID() != INT_NODATA))
578 {
579 bSpanBelow = false;
580 }
581
582 nX++;
583 }
584 }
585
586 // Store this polygon's stored unconsolidated sediment depths
587 pPolygon->SetPreExistingUnconsFine(dStoredUnconsFine);
588 pPolygon->SetPreExistingUnconsSand(dStoredUnconsSand);
589 pPolygon->SetPreExistingUnconsCoarse(dStoredUnconsCoarse);
590
591 // Store this polygon's values for unconsolidated sediment derived from sediment input event(s)
592 pPolygon->SetSedimentInputUnconsFine(dSedimentInputFine);
593 pPolygon->SetSedimentInputUnconsSand(dSedimentInputSand);
594 pPolygon->SetSedimentInputUnconsCoarse(dSedimentInputCoarse);
595
596 // Store this polygon's stored consolidated sediment depths
597 pPolygon->SetPreExistingConsFine(dStoredConsFine);
598 pPolygon->SetPreExistingConsSand(dStoredConsSand);
599 pPolygon->SetPreExistingConsCoarse(dStoredConsCoarse);
600
601 // Store the number of cells in the interior of the polygon
602 pPolygon->SetNumCellsInPolygon(nCellsInPolygon);
603 // LogStream << m_ulIter << ": in MarkPolygonCells() N cells = " << nCellsInPolygon << " in polygon " << nPolyID << endl;
604
605 // Calculate the total volume of seawater on the polygon (m3) and store it
606 double const dSeaVolume = dTotDepth * m_dCellSide;
607 pPolygon->SetSeawaterVolume(dSeaVolume);
608 }
609 }
610
611 // // DEBUG CODE ===========================================================================================================
612 // if (m_ulIter == 109)
613 // {
614 // string strOutFile = m_strOutPath + "00_polygon_fill_";
615 // strOutFile += to_string(m_ulIter);
616 // strOutFile += ".tif";
617 //
618 // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
619 // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
620 // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
621 // pDataSet->SetGeoTransform(m_dGeoTransform);
622 // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
623 // int n = 0;
624 //
625 // int nNumPoly = m_VCoast[0].nGetNumPolygons();
626 // vector<int> nVPerPoly(nNumPoly, 0);
627 // int nNotInPoly = 0;
628 //
629 // for (int nY = 0; nY < m_nYGridSize; nY++)
630 // {
631 // for (int nX = 0; nX < m_nXGridSize; nX++)
632 // {
633 // int nID = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
634 // if (nID == INT_NODATA)
635 // nNotInPoly++;
636 // else
637 // nVPerPoly[nID]++;
638 //
639 // pdRaster[n++] = nID;
640 // }
641 // }
642 //
643 // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
644 // pBand->SetNoDataValue(m_dMissingValue);
645 // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
646 // if (nRet == CE_Failure)
647 // return;
648 //
649 // GDALClose(pDataSet);
650 // delete[] pdRaster;
651 //
652 // for (int nn = 0; nn < nNumPoly; nn++)
653 // LogStream << m_ulIter << ": polygon " << nn << " has " << nVPerPoly[nn] << " cells" << endl;
654 // LogStream << m_ulIter << " Number of cells not in any polygon = " << nNotInPoly << endl;
655 // LogStream << "==================" << endl;
656 // }
657 // // DEBUG CODE ===========================================================================================================
658}
659
660//===============================================================================================================================
662//===============================================================================================================================
664{
665 // Do this for every coast
666 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
667 {
668 int const nNumPolygons = m_VCoast[nCoast].nGetNumPolygons();
669
670 // // DEBUG CODE =================
671 // for (int m = 0; m < m_VCoast[nCoast].nGetNumProfiles(); m++)
672 // {
673 // CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(m);
674 //
675 // LogStream << m << "\t" << pProfile->nGetProfileID() << "\t";
676 //
677 // int nPointsInProfile = pProfile->nGetProfileSize();
678 //
679 // for (int nPoint = 0; nPoint < nPointsInProfile; nPoint++)
680 // {
681 // CGeom2DPoint Pt = *pProfile->pPtGetPointInProfile(nPoint);
682 // LogStream << "{" << Pt.dGetX() << ", " << Pt.dGetY() << "}";
683 // }
684 // LogStream << endl;
685 // }
686 // // DEBUG CODE =================
687
688 // Do this for every coastal polygon, in along-coast sequence
689 for (int nn = 0; nn < nNumPolygons; nn++)
690 {
691 CGeomCoastPolygon* pThisPolygon = m_VCoast[nCoast].pGetPolygon(nn);
692 int const nThisPolygon = pThisPolygon->nGetPolygonCoastID();
693
694 vector<int> nVUpCoastAdjacentPolygon;
695 vector<int> nVDownCoastAdjacentPolygon;
696
697 vector<double> dVUpCoastBoundaryShare;
698 vector<double> dVDownCoastBoundaryShare;
699
700 // First deal with down-coast adjacent polygons
701 if (pThisPolygon->bIsCoastEndPolygon())
702 {
703 // We are at the end of the coastline, no other polygon is adjacent to the down-coast profile of the end-of-coast polygon
704 nVDownCoastAdjacentPolygon.push_back(INT_NODATA);
705 dVDownCoastBoundaryShare.push_back(1);
706
707 // Store in the polygon
708 pThisPolygon->SetDownCoastAdjacentPolygons(&nVDownCoastAdjacentPolygon);
709 pThisPolygon->SetDownCoastAdjacentPolygonBoundaryShares(&dVDownCoastBoundaryShare);
710 }
711 else
712 {
713 // We are not at the end of the coastline, so there is at least one other polygon adjacent to the down-coast profile of this polygon
714 int const nDownCoastProfile = pThisPolygon->nGetDownCoastProfile();
715 int const nPointsInProfile = pThisPolygon->nGetNumPointsUsedDownCoastProfile();
716
717 double dDownCoastTotBoundaryLen = 0;
718
719 CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(nDownCoastProfile);
720
721 for (int nPoint = 0; nPoint < nPointsInProfile - 1; nPoint++)
722 {
723 CGeom2DPoint const PtStart = *pProfile->pPtGetPointInProfile(nPoint);
724 CGeom2DPoint const PtEnd = *pProfile->pPtGetPointInProfile(nPoint + 1);
725
726 // Calculate the length of this segment of the normal profile. Note that it should not be zero, since we checked for duplicate points when creating profiles
727 double const dDistBetween = dGetDistanceBetween(&PtStart, &PtEnd);
728
729 // Find out which polygon is adjacent to each line segment of the polygon's down-coast profile boundary. The basic approach used is to count the number of coincident profiles in each line segment, and (because we are going down-coast) add this number to 'this' polygon's number. However, some of these coincident profiles may be invalid, so we must count only the valid co-incident profiles
730 int const nNumCoinc = pProfile->nGetNumCoincidentProfilesInLineSegment(nPoint);
731
732 // Safety check
733 if (nNumCoinc < 0)
734 continue;
735
736 int nNumValidCoinc = 0;
737
738 for (int nCoinc = 0; nCoinc < nNumCoinc; nCoinc++)
739 {
740 int const nProf = pProfile->nGetCoincidentProfileForLineSegment(nPoint, nCoinc);
741
742 // Safety check
743 if (nProf == -1)
744 continue;
745
746 CGeomProfile const* pProf = m_VCoast[nCoast].pGetProfile(nProf);
747
748 if (pProf->bProfileOK())
749 nNumValidCoinc++;
750 }
751
752 // First stab at calculating the number of the adjacent polygon
753 int nAdj = nThisPolygon + nNumValidCoinc;
754
755 // However, if 'this' polygon is close to the end of the coastline, we get polygon numbers greater than the number of polygons i.e. beyond the end of the coastline. If this happens, set the adjacent polygon to 'off-edge'
756 if (nAdj >= nNumPolygons)
757 nAdj = INT_NODATA;
758
759 // Safety check: is this adjacent polygon already present in the list of down-coast adjacent polygons?
760 if (pThisPolygon->bDownCoastIsAlreadyPresent(nAdj))
761 continue;
762
763 // Not already present
764 nVDownCoastAdjacentPolygon.push_back(nAdj);
765
766 dDownCoastTotBoundaryLen += dDistBetween;
767 dVDownCoastBoundaryShare.push_back(dDistBetween);
768 }
769
770 // Calculate the down-coast boundary share
771 for (unsigned int n = 0; n < dVDownCoastBoundaryShare.size(); n++)
772 {
773 // Safety check
774 if (bFPIsEqual(dDownCoastTotBoundaryLen, 0.0, TOLERANCE))
775 {
776 nVDownCoastAdjacentPolygon.pop_back();
777 dVDownCoastBoundaryShare.pop_back();
778 continue;
779 }
780
781 dVDownCoastBoundaryShare[n] /= dDownCoastTotBoundaryLen;
782 }
783
784 // Store in the polygon
785 pThisPolygon->SetDownCoastAdjacentPolygons(&nVDownCoastAdjacentPolygon);
786 pThisPolygon->SetDownCoastAdjacentPolygonBoundaryShares(&dVDownCoastBoundaryShare);
787 }
788
789 // Now deal with up-coast adjacent polygons
790 if (pThisPolygon->bIsCoastStartPolygon())
791 {
792 // We are at the start of the coastline, no other polygon is adjacent to the up-coast profile of the start-of-coast polygon
793 nVUpCoastAdjacentPolygon.push_back(INT_NODATA);
794 dVUpCoastBoundaryShare.push_back(1);
795
796 // Store in the polygon
797 pThisPolygon->SetUpCoastAdjacentPolygons(&nVUpCoastAdjacentPolygon);
798 pThisPolygon->SetUpCoastAdjacentPolygonBoundaryShares(&dVUpCoastBoundaryShare);
799 }
800
801 else
802 {
803 // We are not at the start of the coastline, so there is at least one other polygon adjacent to the up-coast profile of this polygon
804 int const nUpCoastProfile = pThisPolygon->nGetUpCoastProfile();
805 int const nPointsInProfile = pThisPolygon->nGetNumPointsUsedUpCoastProfile();
806
807 double dUpCoastTotBoundaryLen = 0;
808
809 CGeomProfile* pProfile = m_VCoast[nCoast].pGetProfile(nUpCoastProfile);
810
811 for (int nPoint = 0; nPoint < nPointsInProfile - 1; nPoint++)
812 {
813 CGeom2DPoint const PtStart = *pProfile->pPtGetPointInProfile(nPoint);
814 CGeom2DPoint const PtEnd = *pProfile->pPtGetPointInProfile(nPoint + 1);
815
816 // Safety check
817 if (PtStart == PtEnd)
818 // Should never get here, since we checked for duplicate points when creating profiles
819 continue;
820
821 // Calculate the length of this segment of the normal profile
822 double const dDistBetween = dGetDistanceBetween(&PtStart, &PtEnd);
823
824 // Find out which polygon is adjacent to each line segment of the polygon's up-coast profile boundary. The basic approach used is to count the number of coincident profiles in each line segment, and (because we are going up-coast) subtract this number from 'this' polygon's number. However, some of these coincident profiles may be invalid, so we must count only the valid co-incident profiles
825 int const nNumCoinc = pProfile->nGetNumCoincidentProfilesInLineSegment(nPoint);
826
827 // Safety check
828 if (nNumCoinc < 0)
829 continue;
830
831 int nNumValidCoinc = 0;
832
833 for (int nCoinc = 0; nCoinc < nNumCoinc; nCoinc++)
834 {
835 int const nProf = pProfile->nGetCoincidentProfileForLineSegment(nPoint, nCoinc);
836
837 // Safety check
838 if (nProf == -1)
839 continue;
840
841 CGeomProfile const* pProf = m_VCoast[nCoast].pGetProfile(nProf);
842
843 if (pProf->bProfileOK())
844 nNumValidCoinc++;
845 }
846
847 // First stab at calculating the number of the adjacent polygon
848 int nAdj = nThisPolygon - nNumValidCoinc;
849
850 // However, if 'this' polygon is close to the start of the coastline, we get polygon numbers below zero i.e. beyond the start of the coastline. If this happens, set the adjacent polygon to 'off-edge'
851 if (nAdj < 0)
852 nAdj = INT_NODATA;
853
854 // Safety check: is this adjacent polygon already present in the list of up-coast adjacent polygons?
855 if (pThisPolygon->bUpCoastIsAlreadyPresent(nAdj))
856 continue;
857
858 // Not already present
859 nVUpCoastAdjacentPolygon.push_back(nAdj);
860
861 dUpCoastTotBoundaryLen += dDistBetween;
862 dVUpCoastBoundaryShare.push_back(dDistBetween);
863 }
864
865 // Calculate the up-coast boundary share
866 for (unsigned int n = 0; n < dVUpCoastBoundaryShare.size(); n++)
867 {
868 // Safety check
869 if (bFPIsEqual(dUpCoastTotBoundaryLen, 0.0, TOLERANCE))
870 {
871 nVUpCoastAdjacentPolygon.pop_back();
872 dVUpCoastBoundaryShare.pop_back();
873 continue;
874 }
875
876 dVUpCoastBoundaryShare[n] /= dUpCoastTotBoundaryLen;
877 }
878
879 // Store in the polygon
880 pThisPolygon->SetUpCoastAdjacentPolygons(&nVUpCoastAdjacentPolygon);
881 pThisPolygon->SetUpCoastAdjacentPolygonBoundaryShares(&dVUpCoastBoundaryShare);
882 }
883
884 // Finally, calculate the distance between the coast node and the antinode of the polygon
885 double const dPolygonSeawardLen = dGetDistanceBetween(pThisPolygon->pPtiGetNode(), pThisPolygon->pPtiGetAntiNode());
886
887 // And store it
888 pThisPolygon->SetLength(dPolygonSeawardLen);
889
890 // // DEBUG CODE ======================================================================================================================
891 // assert(dVUpCoastBoundaryShare.size() == nVUpCoastAdjacentPolygon.size());
892 // assert(dVDownCoastBoundaryShare.size() == nVDownCoastAdjacentPolygon.size());
893 //
894 // LogStream << m_ulIter << ": polygon = " << nPoly << (pPolygon->bIsPointed() ? " IS TRIANGULAR" : "") << endl;
895 // LogStream << m_ulIter << ": coast " << nCoast << " polygon " << nThisPolygon << endl;
896 //
897 // int nUpCoastProfile = pThisPolygon->nGetUpCoastProfile();
898 // CGeomProfile* pUpCoastProfile = m_VCoast[0].pGetProfile(nUpCoastProfile);
899 // int nUpCoastProfileCells = pUpCoastProfile->nGetNumCellsInProfile();
900 //
901 // int nDownCoastProfile = pThisPolygon->nGetDownCoastProfile();
902 // CGeomProfile* pDownCoastProfile = m_VCoast[0].pGetProfile(nDownCoastProfile);
903 // int nDownCoastProfileCells = pDownCoastProfile->nGetNumCellsInProfile();
904 //
905 // LogStream << "\tUp-coast profile = " << nUpCoastProfile << " down-coast profile = " << nDownCoastProfile << endl;
906 // LogStream << "\tN cells in up-coast profile = " << nUpCoastProfileCells << " N cells in down-coast profile = " << nDownCoastProfileCells << endl;
907 //
908 // LogStream << "\tThere are " << nVUpCoastAdjacentPolygon.size() << " UP-COAST adjacent polygon(s) = ";
909 // for (unsigned int n = 0; n < nVUpCoastAdjacentPolygon.size(); n++)
910 // LogStream << nVUpCoastAdjacentPolygon[n] << " ";
911 // LogStream << endl;
912 //
913 // LogStream << "\tThere are " << nVDownCoastAdjacentPolygon.size() << " DOWN-COAST adjacent polygon(s) = ";
914 // for (unsigned int n = 0; n < nVDownCoastAdjacentPolygon.size(); n++)
915 // LogStream << nVDownCoastAdjacentPolygon[n] << " ";
916 // LogStream << endl;
917 //
918 // double dUpCoastTotBoundaryLen = 0;
919 // LogStream << "\tUP-COAST boundary share(s) = ";
920 // for (unsigned int n = 0; n < dVUpCoastBoundaryShare.size(); n++)
921 // {
922 // dUpCoastTotBoundaryLen += dVUpCoastBoundaryShare[n];
923 // LogStream << dVUpCoastBoundaryShare[n] << " ";
924 // }
925 // LogStream << endl;
926 // LogStream << "\tTotal UP-COAST boundary length = " << dUpCoastTotBoundaryLen << endl;
927 //
928 // double dDownCoastTotBoundaryLen = 0;
929 // LogStream << "\tDOWN-COAST boundary share(s) = ";
930 // for (unsigned int n = 0; n < dVDownCoastBoundaryShare.size(); n++)
931 // {
932 // dDownCoastTotBoundaryLen += dVDownCoastBoundaryShare[n];
933 // LogStream << dVDownCoastBoundaryShare[n] << " ";
934 // }
935 // LogStream << endl;
936 // LogStream << "\tTotal DOWN-COAST boundary length = " << dDownCoastTotBoundaryLen << endl;
937 // // DEBUG CODE ======================================================================================================================
938 }
939 }
940
941 // // DEBUG CODE =================================================================================
942 // for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
943 // {
944 // int nNumPoly = m_VCoast[nCoast].nGetNumPolygons();
945 // for (int nPoly = 0; nPoly < nNumPoly; nPoly++)
946 // {
947 // CGeomCoastPolygon const* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
948 //
949 // int nNumUpCoastPolygons = pPolygon->nGetNumUpCoastAdjacentPolygons();
950 // LogStream << "Polygon " << nPoly << " up-coast polygon(s) = ";
951 // for (int nAdj = 0; nAdj < nNumUpCoastPolygons; nAdj++)
952 // {
953 // int nAdjPoly = pPolygon->nGetUpCoastAdjacentPolygon(nAdj);
954 // LogStream << nAdjPoly;
955 // if (nAdjPoly > nNumPoly-1)
956 // LogStream << "***";
957 // }
958 //
959 // int nNumDownCoastPolygons = pPolygon->nGetNumDownCoastAdjacentPolygons();
960 // LogStream << " down-coast polygon(s) = ";
961 // for (int nAdj = 0; nAdj < nNumDownCoastPolygons; nAdj++)
962 // {
963 // int nAdjPoly = pPolygon->nGetDownCoastAdjacentPolygon(nAdj);
964 // LogStream << nAdjPoly;
965 // if (nAdjPoly > nNumPoly-1)
966 // LogStream << "***";
967 // }
968 // LogStream << endl;
969 // }
970 // }
971 // // DEBUG CODE =================================================================================
972
973 return RTN_OK;
974}
975
976//===============================================================================================================================
978//===============================================================================================================================
979bool CSimulation::bIsWithinPolygon(CGeom2DPoint const* pPtStart, vector<CGeom2DPoint> const* pPtPoints)
980{
981 bool bOddNodes = false;
982
983 int const nPolyCorners = static_cast<int>(pPtPoints->size());
984 int j = nPolyCorners - 1;
985
986 double const dX = pPtStart->dGetX();
987 double const dY = pPtStart->dGetY();
988
989 for (int i = 0; i < nPolyCorners; i++)
990 {
991 double const dCorneriX = pPtPoints->at(i).dGetX();
992 double const dCorneriY = pPtPoints->at(i).dGetY();
993 double const dCornerjX = pPtPoints->at(j).dGetX();
994 double const dCornerjY = pPtPoints->at(j).dGetY();
995
996 if ((dCorneriY < dY && dCornerjY >= dY) || (dCornerjY < dY && dCorneriY >= dY))
997 {
998 if (dCorneriX + (dY - dCorneriY) / (dCornerjY - dCorneriY) * (dCornerjX - dCorneriX) < dX)
999 {
1000 bOddNodes = ! bOddNodes;
1001 }
1002 }
1003
1004 j = i;
1005 }
1006
1007 return bOddNodes;
1008}
1009
1010//===============================================================================================================================
1012//===============================================================================================================================
1013CGeom2DPoint CSimulation::PtFindPointInPolygon(vector<CGeom2DPoint> const* pPtPoints, int const nStartPoint)
1014{
1015 int const nPolySize = static_cast<int>(pPtPoints->size());
1016 int nOffSet = 0;
1017 CGeom2DPoint PtStart;
1018
1019 do
1020 {
1021 // Choose three consecutive points from the polygon
1022 vector<CGeom2DPoint> nVTestPoints;
1023
1024 for (int n = 0; n < 3; n++)
1025 {
1026 int nIndex = n + nStartPoint + nOffSet;
1027
1028 if (nIndex > nPolySize - 1)
1029 nIndex -= nPolySize;
1030
1031 // Safety check
1032 if (nIndex < 0)
1034
1035 nVTestPoints.push_back(pPtPoints->at(nIndex));
1036 }
1037
1038 // Increment ready for next time
1039 nOffSet++;
1040
1041 // Safety check
1042 if (nOffSet >= (nPolySize + 3))
1044
1045 // Check if the halfway point between the first and the third point is inside the polygon
1046 PtStart = PtAverage(&nVTestPoints[0], &nVTestPoints[2]);
1047 } while (! bIsWithinPolygon(&PtStart, pPtPoints));
1048
1049 return PtStart;
1050}
Contains CGeom2DPoint definitions.
Contains CGeom2DIPoint definitions.
Geometry class used to represent 2D point objects with integer coordinates.
Definition 2di_point.h:29
int nGetY(void) const
Returns the CGeom2DIPoint object's integer Y coordinate.
Definition 2di_point.cpp:55
int nGetX(void) const
Returns the CGeom2DIPoint object's integer X coordinate.
Definition 2di_point.cpp:49
Geometry class used to represent 2D point objects with floating-point coordinates.
Definition 2d_point.h:27
double dGetY(void) const
Returns the CGeom2DPoint object's double Y coordinate.
Definition 2d_point.cpp:53
double dGetX(void) const
Returns the CGeom2DPoint object's double X coordinate.
Definition 2d_point.cpp:47
Geometry class used for coast polygon objects.
bool bIsCoastStartPolygon(void) const
Is this polygon the coast-start polygon?
int nGetNumPointsUsedDownCoastProfile(void) const
Return the number of points in the down-coast profile.
CGeom2DIPoint * pPtiGetAntiNode(void)
Get the anti-node (raster grid CRS) which is at other (seaward) end of the polygon from the node.
bool bIsCoastEndPolygon(void) const
Is this polygon the coast-end polygon?
void SetSedimentInputUnconsSand(double const)
Set the value of sand sediment on the polygon derived from sediment input events(s)
void SetPreExistingUnconsFine(double const)
Set the value of pre-existing unconsolidated fine sediment stored on this polygon.
void SetSedimentInputUnconsFine(double const)
Set the value of fine sediment on the polygon derived from sediment input events(s)
void SetPreExistingConsFine(double const)
Set the value of pre-existing consolidated fine sediment stored on this polygon.
void SetPreExistingUnconsSand(double const)
Set the value of pre-existing unconsolidated sand sediment stored on this polygon.
void SetSeawaterVolume(const double)
Set the volume of seawater in the coast polygon.
bool bUpCoastIsAlreadyPresent(int const)
Searches the list of up-coast adjacent polygons for a given polygon number.
int nGetNumPointsUsedUpCoastProfile(void) const
Return the number of points in the up-coast profile.
CGeom2DIPoint PtiGetFillStartPoint(void)
void SetSedimentInputUnconsCoarse(double const)
Set the value of coarse sediment on the polygon derived from sediment input events(s)
void SetDownCoastAdjacentPolygonBoundaryShares(vector< double > const *)
Sets the boundary shares for all down-coast adjacent polygons.
int nGetPolygonCoastID(void) const
Get the coast ID, this is the same as the down-coast sequence of polygons.
int nGetUpCoastProfile(void) const
Return the number of the up-coast profile.
void SetLength(double const)
Sets the polygon's length.
void SetPreExistingConsCoarse(double const)
Set the value of pre-existing consolidated coarse sediment stored on this polygon.
void SetPreExistingUnconsCoarse(double const)
Set the value of pre-existing unconsolidated coarse sediment stored on this polygon.
void AppendVertex(CGeom2DIPoint const *)
Appends the point cordinates (grid CRS) for a polygon vertex.
void SetDownCoastAdjacentPolygons(vector< int > const *)
Sets all down-coast adjacent polygons.
void SetUpCoastAdjacentPolygonBoundaryShares(vector< double > const *)
Sets the boundary shares for all up-coast adjacent polygons.
void SetUpCoastAdjacentPolygons(vector< int > const *)
Sets all up-coast adjacent polygons.
void SetPreExistingConsSand(double const)
Set the value of pre-existing consolidated sand sediment stored on this polygon.
int nGetDownCoastProfile(void) const
Return the number of the down-coast profile.
bool bDownCoastIsAlreadyPresent(int const)
Searches the list of down-coast adjacent polygons for a given polygon number.
void SetNumCellsInPolygon(int const)
Set the number of cells in the polygon.
CGeom2DIPoint * pPtiGetNode(void)
Get the grid coordinates of the cell on which the node sits.
int nGetNumCoincidentProfilesInLineSegment(int const)
Returns the count of coincident profiles in a specified line segment, or -1 if the line segment does ...
bool bFindProfileInCoincidentProfiles(int const)
Returns true if the given profile number is one of the coincident profiles of the a specified line se...
int nGetCoincidentProfileForLineSegment(int const, int const) const
Returns the numbers of coincident profiles.
void GetMostCoastwardSharedLineSegment(int const, int &, int &)
Finds the number of the most coastward line segment for which the two profiles are coincident,...
Geometry class used to represent coast profile objects.
Definition profile.h:37
int nGetProfileID(void) const
Returns the profile's this-coast ID.
Definition profile.cpp:80
int nGetCoastPoint(void) const
Returns the coast point at which the profile starts.
Definition profile.cpp:86
CGeomProfile * pGetDownCoastAdjacentProfile(void) const
Returns the down-coast adjacent profile, returns NULL if there is no down-coast adjacent profile.
Definition profile.cpp:491
CGeom2DIPoint * pPtiGetCellInProfile(int const)
Returns a single cell (grid CRS) in the profile.
Definition profile.cpp:521
int nGetNumCellsInProfile(void) const
Returns the number of cells in the profile.
Definition profile.cpp:535
bool bProfileOK(void) const
Returns true if this is a problem-free profile, and is not a start-of-coast and is not an end-of-coas...
Definition profile.cpp:230
CGeom2DIPoint * pPtiGetEndPoint(void)
Returns a pointer to the location of the cell (grid CRS) on which the profile ends.
Definition profile.cpp:104
CGeom2DPoint * pPtGetPointInProfile(int const)
Returns a single point (external CRS) from the profile.
Definition profile.cpp:367
bool bOKIncStartAndEndOfCoast(void) const
Returns true if this is a problem-free profile (however it could be a start-of-coast or an end-of-coa...
Definition profile.cpp:265
int nGetProfileSize(void) const
Returns the number of external CRS points in the profile (only two, initally; and always just two for...
Definition profile.cpp:361
CGeom2DIPoint * pPtiGetStartPoint(void)
Returns a pointer to the location of the cell (grid CRS) on which the profile starts.
Definition profile.cpp:92
bool bStartOfCoast(void) const
Returns the switch to indicate whether this is a start-of-coast profile.
Definition profile.cpp:116
bool bEndOfCoast(void) const
Returns the switch to indicate whether this is an end-of-coast profile.
Definition profile.cpp:128
CGeomRasterGrid * m_pRasterGrid
Pointer to the raster grid object.
int m_nXGridSize
The size of the grid in the x direction.
Definition simulation.h:457
CGeom2DIPoint PtiExtCRSToGridRound(CGeom2DPoint const *) const
Transforms a pointer to a CGeom2DPoint in the external CRS to the equivalent CGeom2DIPoint in the ras...
ofstream LogStream
vector< CRWCoast > m_VCoast
The coastline objects.
int m_nYGridSize
The size of the grid in the y direction.
Definition simulation.h:460
static double dGetDistanceBetween(CGeom2DPoint const *, CGeom2DPoint const *)
Returns the distance (in external CRS) between two points.
static CGeom2DPoint PtAverage(CGeom2DPoint const *, CGeom2DPoint const *)
Returns a point (external CRS) which is the average of (i.e. is midway between) two other external CR...
double m_dStartIterSuspFineInPolygons
Depth (m) of fine suspended sediment at the start of the simulation (only cells in polygons)
Definition simulation.h:985
void KeepWithinValidGrid(int &, int &) const
int nCreateAllPolygons(void)
Create polygons, and mark the polygon boundaries on the raster grid.
void RasterizePolygonJoiningLine(int const, CGeom2DIPoint const *, CGeom2DIPoint const *, int const)
Puts a polygon 'joining line' (the line which is the seaward boundary of the polygon,...
static CGeom2DPoint PtFindPointInPolygon(vector< CGeom2DPoint > const *, int const)
Finds a point in a polygon: is guaranteed to succeed, as every strictly closed polygon has at least o...
void MarkPolygonCells(void)
Marks cells of the raster grid that are within each coastal polygon. The cell-by-cell fill (aka 'floo...
bool bIsWithinValidGrid(int const, int const) const
unsigned long m_ulIter
The number of the current iteration (time step)
Definition simulation.h:598
int m_nNumPolygonGlobal
Number of global (all coasts) polygons.
Definition simulation.h:523
double m_dCellSide
Length of a cell side (in external CRS units)
Definition simulation.h:658
int nDoPolygonSharedBoundaries(void)
For between-polygon potential sediment routing: find which are the adjacent polygons,...
static bool bIsWithinPolygon(CGeom2DPoint const *, vector< CGeom2DPoint > const *)
Determines whether a point is within a polygon: however if the point is exactly on the edge of the po...
vector< CGeomCoastPolygon * > m_pVCoastPolygon
Pointers to coast polygon objects, in down-coast sequence TODO 044 Will need to use global polygon ID...
This file contains global definitions for CoastalME.
int const NO_NONZERO_THICKNESS_LAYERS
Definition cme.h:778
int const INT_NODATA
Definition cme.h:476
T tMin(T a, T b)
Definition cme.h:1265
double const TOLERANCE
Definition cme.h:823
string const ERR
Definition cme.h:902
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1303
T tMax(T a, T b)
Definition cme.h:1252
int const RTN_ERR_BAD_MULTILINE
Definition cme.h:740
int const RTN_OK
Definition cme.h:694
double const DBL_NODATA
Definition cme.h:834
T tAbs(T a)
Definition cme.h:1277
Contains CRWCoast definitions.
Contains CSimulation definitions.
int nRound(double const d)
Version of the above that returns an int.