gis_raster.cpp

Go to the documentation of this file.

 
/* ===============================================================================================================================
 
   This file is part of CoastalME, the Coastal Modelling Environment.
 
   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.
 
   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.
 
   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.
 
===============================================================================================================================*/
#include <assert.h>
 
#include <cstdio>
 
#include <cmath>
using std::hypot;
using std::isnan;
using std::sqrt;
using std::atan2;
using std::isfinite;
 
#include <vector>
using std::vector;
 
#include <iostream>
using std::cerr;
// using std::cout;
using std::endl;
using std::ios;
 
#include <fstream>
using std::ifstream;
 
#include <sstream>
using std::stringstream;
 
#include <string>
using std::to_string;
 
#include <gdal.h>
#include <gdal_priv.h>
#include <gdal_alg.h>
#include <cpl_conv.h>
#include <cpl_error.h>
#include <cpl_string.h>
 
#include "cme.h"
#include "simulation.h"
#include "coast.h"
#include "2di_point.h"
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::InitializeGDALPerformance(void)
{
   // Configure GDAL for optimal performance
   // Enable GDAL threading - use all available CPU cores
#ifdef _OPENMP
   CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
#else
   CPLSetConfigOption("GDAL_NUM_THREADS", "2"); // Fallback for non-OpenMP builds
#endif
 
   // Optimize GDAL memory usage and caching
   CPLSetConfigOption("GDAL_CACHEMAX", "1024");                // 1GB cache for large grids
   CPLSetConfigOption("GDAL_DISABLE_READDIR_ON_OPEN", "TRUE"); // Faster file access
   CPLSetConfigOption("VSI_CACHE", "TRUE");                    // Enable virtual file system cache
   CPLSetConfigOption("VSI_CACHE_SIZE", "256000000");          // 256MB VSI cache
 
   // Optimize grid creation performance
   CPLSetConfigOption("GDAL_GRID_MAX_POINTS_PER_QUADTREE_LEAF", "512"); // Faster spatial indexing
 
   // Disable GDAL warnings for cleaner output (optional)
   // CPLSetConfigOption("CPL_LOG", "/dev/null");
 
   LogStream << "GDAL performance optimizations enabled" << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nReadRasterBasementDEM(void)
{
   // Initialize GDAL performance settings (only needs to be done once)
   static bool bGDALInitialized = false;
 
   if (! bGDALInitialized)
   {
      InitializeGDALPerformance();
      bGDALInitialized = true;
   }
 
   // Use GDAL to create a dataset object, which then opens the DEM file
   GDALDataset *pGDALDataset = static_cast<GDALDataset*>(GDALOpen(m_strInitialBasementDEMFile.c_str(), GA_ReadOnly));
 
   if (NULL == pGDALDataset)
   {
      // Can't open file (note will already have sent GDAL error message to stdout)
      cerr << ERR << "cannot open " << m_strInitialBasementDEMFile << " for input: " << CPLGetLastErrorMsg() << endl;
      return RTN_ERR_DEMFILE;
   }
 
   // Opened OK, so get GDAL basement DEM dataset information
   m_strGDALBasementDEMDriverCode = pGDALDataset->GetDriver()->GetDescription();
   m_strGDALBasementDEMDriverDesc = pGDALDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME);
   m_strGDALBasementDEMProjection = pGDALDataset->GetProjectionRef();
 
   // If we have reference units, then check that they are in metres
   if (! m_strGDALBasementDEMProjection.empty())
   {
      string const strTmp = strToLower(&m_strGDALBasementDEMProjection);
 
      if ((strTmp.find("meter") == string::npos) && (strTmp.find("metre") == string::npos))
      {
         // error: x-y values must be in metres
         cerr << ERR << "GIS file x-y values (" << m_strGDALBasementDEMProjection << ") in " << m_strInitialBasementDEMFile << " must be in metres" << endl;
         return RTN_ERR_DEMFILE;
      }
   }
 
   // Now get dataset size, and do some rudimentary checks
   m_nXGridSize = pGDALDataset->GetRasterXSize();
 
   if (m_nXGridSize == 0)
   {
      // Error: silly number of columns specified
      cerr << ERR << "invalid number of columns (" << m_nXGridSize << ") in " << m_strInitialBasementDEMFile << endl;
      return RTN_ERR_DEMFILE;
   }
 
   m_nYGridSize = pGDALDataset->GetRasterYSize();
 
   if (m_nYGridSize == 0)
   {
      // Error: silly number of rows specified
      cerr << ERR << "invalid number of rows (" << m_nYGridSize << ") in " << m_strInitialBasementDEMFile << endl;
      return RTN_ERR_DEMFILE;
   }
 
   // Get geotransformation info (see http://www.gdal.org/classGDALDataset.html)
   if (CE_Failure == pGDALDataset->GetGeoTransform(m_dGeoTransform))
   {
      // Can't get geotransformation (note will already have sent GDAL error message to stdout)
      cerr << ERR << CPLGetLastErrorMsg() << " in " << m_strInitialBasementDEMFile << endl;
      return RTN_ERR_DEMFILE;
   }
 
   // CoastalME can only handle rasters that are oriented N-S and W-E. (If you need to work with a raster that is oriented differently, then you must rotate it before running CoastalME). So here we check whether row rotation (m_dGeoTransform[2]) and column rotation (m_dGeoTransform[4]) are both zero. See https://gdal.org/tutorials/geotransforms_tut.html
   if ((! bFPIsEqual(m_dGeoTransform[2], 0.0, TOLERANCE)) || (! bFPIsEqual(m_dGeoTransform[4], 0.0, TOLERANCE)))
   {
      // Error: not oriented NS and W-E
      cerr << ERR << m_strInitialBasementDEMFile << " is not oriented N-S and W-E. Row rotation = " << m_dGeoTransform[2] << " and column rotation = " << m_dGeoTransform[4] << endl;
      return (RTN_ERR_RASTER_FILE_READ);
   }
 
   // Get the X and Y cell sizes, in external CRS units. Note that while the cell is supposed to be square, it may not be exactly so due to oddities with some GIS calculations
   double const dCellSideX = tAbs(m_dGeoTransform[1]);
   double const dCellSideY = tAbs(m_dGeoTransform[5]);
 
   // Check that the cell is more or less square
   if (! bFPIsEqual(dCellSideX, dCellSideY, 1e-2))
   {
      // Error: cell is not square enough
      cerr << ERR << "cell is not square in " << m_strInitialBasementDEMFile << ", is " << dCellSideX << " x " << dCellSideY << endl;
      return (RTN_ERR_RASTER_FILE_READ);
   }
 
   // Calculate the average length of cell side, the cell's diagonal, and the area of a cell (in external CRS units)
   m_dCellSide = (dCellSideX + dCellSideY) / 2.0;
   m_dCellArea = m_dCellSide * m_dCellSide;
   m_dCellDiagonal = hypot(m_dCellSide, m_dCellSide);
 
   // And calculate the inverse values
   m_dInvCellSide = 1 / m_dCellSide;
   m_dInvCellDiagonal = 1 / m_dCellDiagonal;
 
   // Save some values in external CRS
   m_dNorthWestXExtCRS = m_dGeoTransform[0] - (m_dGeoTransform[1] / 2);
   m_dNorthWestYExtCRS = m_dGeoTransform[3] - (m_dGeoTransform[5] / 2);
   m_dSouthEastXExtCRS = m_dGeoTransform[0] + (m_nXGridSize * m_dGeoTransform[1]) + (m_dGeoTransform[1] / 2);
   m_dSouthEastYExtCRS = m_dGeoTransform[3] + (m_nYGridSize * m_dGeoTransform[5]) + (m_dGeoTransform[5] / 2);
 
   // And calc the grid area in external CRS units
   m_dExtCRSGridArea = tAbs(m_dNorthWestXExtCRS - m_dSouthEastXExtCRS) * tAbs(m_dNorthWestYExtCRS * m_dSouthEastYExtCRS);
 
   // Now get GDAL raster band information
   GDALRasterBand *pGDALBand = pGDALDataset->GetRasterBand(1);
   int nBlockXSize = 0, nBlockYSize = 0;
   pGDALBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
   m_strGDALBasementDEMDataType = GDALGetDataTypeName(pGDALBand->GetRasterDataType());
 
   // If we have value units, then check them
   string const strUnits = pGDALBand->GetUnitType();
 
   if ((!strUnits.empty()) && (strUnits.find('m') == string::npos))
   {
      // Error: value units must be m
      cerr << ERR << "DEM vertical units are (" << strUnits << " ) in " << m_strInitialBasementDEMFile << ", should be 'm'" << endl;
      return RTN_ERR_DEMFILE;
   }
 
   // If present, get the missing value (NODATA) setting
   CPLPushErrorHandler(CPLQuietErrorHandler);                // Needed to get next line to fail silently, if it fails
   double const dMissingValue = pGDALBand->GetNoDataValue(); // Will fail for some formats
   CPLPopErrorHandler();
 
   if (! bFPIsEqual(dMissingValue, m_dMissingValue, TOLERANCE))
   {
      cerr << "   " << NOTE << "NODATA value in " << m_strInitialBasementDEMFile << " is " << dMissingValue << "\n         instead using CoastalME's default floating-point NODATA value " << m_dMissingValue << endl;
   }
 
   // Next allocate memory for a 2D array of raster cell objects: tell the user what is happening
   AnnounceAllocateMemory();
   int const nRet = m_pRasterGrid->nCreateGrid();
 
   if (nRet != RTN_OK)
      return nRet;
 
   // Allocate memory for a 1D floating-point array, to hold the scan line for GDAL
   double *pdScanline = new double[m_nXGridSize];
 
   if (NULL == pdScanline)
   {
      // Error, can't allocate memory
      cerr << ERR << "cannot allocate memory for " << m_nXGridSize << " x 1D array" << endl;
      return (RTN_ERR_MEMALLOC);
   }
 
   // Now read in the data
   for (int j = 0; j < m_nYGridSize; j++)
   {
      // Read scanline
      if (CE_Failure == pGDALBand->RasterIO(GF_Read, 0, j, m_nXGridSize, 1, pdScanline, m_nXGridSize, 1, GDT_Float64, 0, 0, NULL))
      {
         // Error while reading scanline
         cerr << ERR << CPLGetLastErrorMsg() << " in " << m_strInitialBasementDEMFile << endl;
         return RTN_ERR_DEMFILE;
      }
 
      // All OK, so read scanline into cell elevations (including any missing values)
      for (int i = 0; i < m_nXGridSize; i++)
      {
         double dTmp = pdScanline[i];
 
         if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
            dTmp = m_dMissingValue;
 
         m_pRasterGrid->m_Cell[i][j].SetBasementElev(dTmp);
      }
   }
 
   // Finished, so get rid of dataset object
   GDALClose(pGDALDataset);
 
   // Get rid of memory allocated to this array
   delete[] pdScanline;
 
   return RTN_OK;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nMarkBoundingBoxEdgeCells(void)
{
   // The bounding box must touch the edge of the grid at least once on each side of the grid, so store these points. Search in a clockwise direction around the edge of the grid
   vector<CGeom2DIPoint> VPtiBoundingBoxCorner;
 
   // Start with the top (north) edge
   bool bFound = false;
 
   for (int nX = 0; nX < m_nXGridSize; nX++)
   {
      if (bFound)
         break;
 
      for (int nY = 0; nY < m_nYGridSize; nY++)
      {
         if (! m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            CGeom2DIPoint const PtiTmp(nX, nY);
            VPtiBoundingBoxCorner.push_back(PtiTmp);
            bFound = true;
            break;
         }
      }
   }
 
   if (! bFound)
   {
      if (m_nLogFileDetail >= LOG_FILE_ALL)
         LogStream << m_ulIter << ": north (top) edge of bounding box not found" << endl;
 
      return RTN_ERR_BOUNDING_BOX;
   }
 
   // Do the same for the right (east) edge
   bFound = false;
 
   for (int nY = 0; nY < m_nYGridSize; nY++)
   {
      if (bFound)
         break;
 
      for (int nX = m_nXGridSize - 1; nX >= 0; nX--)
      {
         if (! m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            CGeom2DIPoint const PtiTmp(nX, nY);
            VPtiBoundingBoxCorner.push_back(PtiTmp);
            bFound = true;
            break;
         }
      }
   }
 
   if (! bFound)
   {
      if (m_nLogFileDetail >= LOG_FILE_ALL)
         LogStream << m_ulIter << ": east (right) edge of bounding box not found" << endl;
 
      return RTN_ERR_BOUNDING_BOX;
   }
 
   // Do the same for the south (bottom) edge
   bFound = false;
 
   for (int nX = m_nXGridSize - 1; nX >= 0; nX--)
   {
      if (bFound)
         break;
 
      for (int nY = m_nYGridSize - 1; nY >= 0; nY--)
      {
         if (! m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            CGeom2DIPoint const PtiTmp(nX, nY);
            VPtiBoundingBoxCorner.push_back(PtiTmp);
            bFound = true;
            break;
         }
      }
   }
 
   if (! bFound)
   {
      if (m_nLogFileDetail >= LOG_FILE_ALL)
         LogStream << m_ulIter << ": south (bottom) edge of bounding box not found" << endl;
 
      return RTN_ERR_BOUNDING_BOX;
   }
 
   // And finally repeat for the west (left) edge
   bFound = false;
 
   for (int nY = m_nYGridSize - 1; nY >= 0; nY--)
   {
      if (bFound)
         break;
 
      for (int nX = 0; nX < m_nXGridSize; nX++)
      {
         if (! m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            CGeom2DIPoint const PtiTmp(nX, nY);
            VPtiBoundingBoxCorner.push_back(PtiTmp);
            bFound = true;
            break;
         }
      }
   }
 
   if (! bFound)
   {
      if (m_nLogFileDetail >= LOG_FILE_ALL)
         LogStream << m_ulIter << ": west (left) edge of bounding box not found" << endl;
 
      return RTN_ERR_BOUNDING_BOX;
   }
 
   // OK, so we have a point on each side of the grid, so start at this point and find the edges of the bounding box. Go round in a clockwise direction: top (north) edge first
   for (int nX = VPtiBoundingBoxCorner[0].nGetX(); nX <= VPtiBoundingBoxCorner[1].nGetX(); nX++)
   {
      bFound = false;
 
      for (int nY = VPtiBoundingBoxCorner[0].nGetY(); nY < m_nYGridSize; nY++)
      {
         if (m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            m_ulMissingValueBasementCells++;
            continue;
         }
 
         // Found a bounding box edge cell
         m_pRasterGrid->m_Cell[nX][nY].SetBoundingBoxEdge(NORTH);
 
         m_VEdgeCell.push_back(CGeom2DIPoint(nX, nY));
         m_VEdgeCellEdge.push_back(NORTH);
 
         bFound = true;
         break;
      }
 
      if (! bFound)
      {
         if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
            LogStream << m_ulIter << ": could not find a bounding box edge cell for grid column " << nX << endl;
 
         return RTN_ERR_BOUNDING_BOX;
      }
   }
 
   // Right (east) edge
   for (int nY = VPtiBoundingBoxCorner[1].nGetY(); nY <= VPtiBoundingBoxCorner[2].nGetY(); nY++)
   {
      bFound = false;
 
      for (int nX = VPtiBoundingBoxCorner[1].nGetX(); nX >= 0; nX--)
      {
         if (m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            m_ulMissingValueBasementCells++;
            continue;
         }
 
         // Found a bounding box edge cell
         m_pRasterGrid->m_Cell[nX][nY].SetBoundingBoxEdge(EAST);
 
         m_VEdgeCell.push_back(CGeom2DIPoint(nX, nY));
         m_VEdgeCellEdge.push_back(EAST);
 
         bFound = true;
         break;
      }
 
      if (! bFound)
      {
         if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
            LogStream << m_ulIter << ": could not find a bounding box edge cell for grid row " << nY << endl;
 
         return RTN_ERR_BOUNDING_BOX;
      }
   }
 
   // Bottom (south) edge
   for (int nX = VPtiBoundingBoxCorner[2].nGetX(); nX >= VPtiBoundingBoxCorner[3].nGetX(); nX--)
   {
      bFound = false;
 
      for (int nY = VPtiBoundingBoxCorner[2].nGetY(); nY >= 0; nY--)
      {
         if (m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            m_ulMissingValueBasementCells++;
            continue;
         }
 
         // Found a bounding box edge cell
         m_pRasterGrid->m_Cell[nX][nY].SetBoundingBoxEdge(SOUTH);
 
         m_VEdgeCell.push_back(CGeom2DIPoint(nX, nY));
         m_VEdgeCellEdge.push_back(SOUTH);
 
         bFound = true;
         break;
      }
 
      if (! bFound)
      {
         if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
            LogStream << m_ulIter << ": could not find a bounding box edge cell for grid column " << nX << endl;
 
         return RTN_ERR_BOUNDING_BOX;
      }
   }
 
   // Left (west) edge
   for (int nY = VPtiBoundingBoxCorner[3].nGetY(); nY >= VPtiBoundingBoxCorner[0].nGetY(); nY--)
   {
      for (int nX = VPtiBoundingBoxCorner[3].nGetX(); nX < m_nXGridSize - 1; nX++)
      // for (int nX = VPtiBoundingBoxCorner[3].nGetX(); nX < VPtiBoundingBoxCorner[3].nGetX(); nX++)
      {
         if (m_pRasterGrid->m_Cell[nX][nY].bBasementElevIsMissingValue())
         {
            m_ulMissingValueBasementCells++;
            continue;
         }
 
         // Found a bounding box edge cell
         m_pRasterGrid->m_Cell[nX][nY].SetBoundingBoxEdge(WEST);
 
         m_VEdgeCell.push_back(CGeom2DIPoint(nX, nY));
         m_VEdgeCellEdge.push_back(WEST);
 
         bFound = true;
         break;
      }
 
      if (! bFound)
      {
         if (m_nLogFileDetail >= LOG_FILE_MIDDLE_DETAIL)
            LogStream << m_ulIter << ": could not find a bounding box edge cell for grid row " << nY << endl;
 
         return RTN_ERR_BOUNDING_BOX;
      }
   }
 
   return RTN_OK;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nReadRasterGISFile(int const nDataItem, int const nLayer)
{
   string
       strGISFile,
       strDriverCode,
       strDriverDesc,
       strProjection,
       strDataType;
 
   switch (nDataItem)
   {
   case (LANDFORM_RASTER):
      // Initial Landform Class GIS data
      strGISFile = m_strInitialLandformFile;
      break;
 
   case (INTERVENTION_CLASS_RASTER):
      // Intervention class
      strGISFile = m_strInterventionClassFile;
      break;
 
   case (INTERVENTION_HEIGHT_RASTER):
      // Intervention height
      strGISFile = m_strInterventionHeightFile;
      break;
 
   case (SUSP_SED_RASTER):
      // Initial Suspended Sediment GIS data
      strGISFile = m_strInitialSuspSedimentFile;
      break;
 
   case (FINE_UNCONS_RASTER):
      // Initial Unconsolidated Fine Sediment GIS data
      strGISFile = m_VstrInitialFineUnconsSedimentFile[nLayer];
      break;
 
   case (SAND_UNCONS_RASTER):
      // Initial Unconsolidated Sand Sediment GIS data
      strGISFile = m_VstrInitialSandUnconsSedimentFile[nLayer];
      break;
 
   case (COARSE_UNCONS_RASTER):
      // Initial Unconsolidated Coarse Sediment GIS data
      strGISFile = m_VstrInitialCoarseUnconsSedimentFile[nLayer];
      break;
 
   case (FINE_CONS_RASTER):
      // Initial Consolidated Fine Sediment GIS data
      strGISFile = m_VstrInitialFineConsSedimentFile[nLayer];
      break;
 
   case (SAND_CONS_RASTER):
      // Initial Consolidated Sand Sediment GIS data
      strGISFile = m_VstrInitialSandConsSedimentFile[nLayer];
      break;
 
   case (COARSE_CONS_RASTER):
      // Initial Consolidated Coarse Sediment GIS data
      strGISFile = m_VstrInitialCoarseConsSedimentFile[nLayer];
      break;
   }
 
   // Do we have a filename for this data item? If we don't then just return
   if (!strGISFile.empty())
   {
      // We do have a filename, so use GDAL to create a dataset object, which then opens the GIS file
      GDALDataset *pGDALDataset = static_cast<GDALDataset*>(GDALOpen(strGISFile.c_str(), GA_ReadOnly));
 
      if (NULL == pGDALDataset)
      {
         // Can't open file (note will already have sent GDAL error message to stdout)
         cerr << ERR << "cannot open " << strGISFile << " for input: " << CPLGetLastErrorMsg() << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      // Opened OK, so get dataset information
      strDriverCode = pGDALDataset->GetDriver()->GetDescription();
      strDriverDesc = pGDALDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME);
      strProjection = pGDALDataset->GetProjectionRef();
 
      // Get geotransformation info
      double dGeoTransform[6];
 
      if (CE_Failure == pGDALDataset->GetGeoTransform(dGeoTransform))
      {
         // Can't get geotransformation (note will already have sent GDAL error message to stdout)
         cerr << ERR << CPLGetLastErrorMsg() << " in " << strGISFile << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      // Now get dataset size, and do some checks
      int const nTmpXSize = pGDALDataset->GetRasterXSize();
 
      if (nTmpXSize != m_nXGridSize)
      {
         // Error: incorrect number of columns specified
         cerr << ERR << "different number of columns in " << strGISFile << " (" << nTmpXSize << ") and " << m_strInitialBasementDEMFile << "(" << m_nXGridSize << ")" << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      int const nTmpYSize = pGDALDataset->GetRasterYSize();
 
      if (nTmpYSize != m_nYGridSize)
      {
         // Error: incorrect number of rows specified
         cerr << ERR << "different number of rows in " << strGISFile << " (" << nTmpYSize << ") and " << m_strInitialBasementDEMFile << " (" << m_nYGridSize << ")" << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      double dTmp = m_dGeoTransform[0] - (m_dGeoTransform[1] / 2);
 
      if (! bFPIsEqual(dTmp, m_dNorthWestXExtCRS, TOLERANCE))
      {
         // Error: different min x from DEM file
         cerr << ERR << "different min x values in " << strGISFile << " (" << dTmp << ") and " << m_strInitialBasementDEMFile << " (" << m_dNorthWestXExtCRS << ")" << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      dTmp = m_dGeoTransform[3] - (m_dGeoTransform[5] / 2);
 
      if (! bFPIsEqual(dTmp, m_dNorthWestYExtCRS, TOLERANCE))
      {
         // Error: different min x from DEM file
         cerr << ERR << "different min y values in " << strGISFile << " (" << dTmp << ") and " << m_strInitialBasementDEMFile << " (" << m_dNorthWestYExtCRS << ")" << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      double const dTmpResX = tAbs(dGeoTransform[1]);
 
      if (! bFPIsEqual(dTmpResX, m_dCellSide, 1e-2))
      {
         // Error: different cell size in X direction: note that due to rounding errors in some GIS packages, must expect some discrepancies
         cerr << ERR << "cell size in X direction (" << dTmpResX << ") in " << strGISFile << " differs from cell size in of basement DEM (" << m_dCellSide << ")" << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      double const dTmpResY = tAbs(dGeoTransform[5]);
 
      if (! bFPIsEqual(dTmpResY, m_dCellSide, 1e-2))
      {
         // Error: different cell size in Y direction: note that due to rounding errors in some GIS packages, must expect some discrepancies
         cerr << ERR << "cell size in Y direction (" << dTmpResY << ") in " << strGISFile << " differs from cell size of basement DEM (" << m_dCellSide << ")" << endl;
         return (RTN_ERR_RASTER_FILE_READ);
      }
 
      // Now get GDAL raster band information
      GDALRasterBand *pGDALBand = pGDALDataset->GetRasterBand(1); // TODO 028 Give a message if there are several bands
      int nBlockXSize = 0, nBlockYSize = 0;
      pGDALBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
      strDataType = GDALGetDataTypeName(pGDALBand->GetRasterDataType());
 
      switch (nDataItem)
      {
      case (LANDFORM_RASTER):
         // Initial Landform Class GIS data
         m_strGDALLDriverCode = strDriverCode;
         m_strGDALLDriverDesc = strDriverDesc;
         m_strGDALLProjection = strProjection;
         m_strGDALLDataType = strDataType;
         break;
 
      case (INTERVENTION_CLASS_RASTER):
         // Intervention class
         m_strGDALICDriverCode = strDriverCode;
         m_strGDALICDriverDesc = strDriverDesc;
         m_strGDALICProjection = strProjection;
         m_strGDALICDataType = strDataType;
         break;
 
      case (INTERVENTION_HEIGHT_RASTER):
         // Intervention height
         m_strGDALIHDriverCode = strDriverCode;
         m_strGDALIHDriverDesc = strDriverDesc;
         m_strGDALIHProjection = strProjection;
         m_strGDALIHDataType = strDataType;
         break;
 
      case (SUSP_SED_RASTER):
         // Initial Suspended Sediment GIS data
         m_strGDALISSDriverCode = strDriverCode;
         m_strGDALISSDriverDesc = strDriverDesc;
         m_strGDALISSProjection = strProjection;
         m_strGDALISSDataType = strDataType;
         break;
 
      case (FINE_UNCONS_RASTER):
         // Initial Unconsolidated Fine Sediment GIS data
         m_VstrGDALIUFDriverCode[nLayer] = strDriverCode;
         m_VstrGDALIUFDriverDesc[nLayer] = strDriverDesc;
         m_VstrGDALIUFProjection[nLayer] = strProjection;
         m_VstrGDALIUFDataType[nLayer] = strDataType;
         break;
 
      case (SAND_UNCONS_RASTER):
         // Initial Unconsolidated Sand Sediment GIS data
         m_VstrGDALIUSDriverCode[nLayer] = strDriverCode;
         m_VstrGDALIUSDriverDesc[nLayer] = strDriverDesc;
         m_VstrGDALIUSProjection[nLayer] = strProjection;
         m_VstrGDALIUSDataType[nLayer] = strDataType;
         break;
 
      case (COARSE_UNCONS_RASTER):
         // Initial Unconsolidated Coarse Sediment GIS data
         m_VstrGDALIUCDriverCode[nLayer] = strDriverCode;
         m_VstrGDALIUCDriverDesc[nLayer] = strDriverDesc;
         m_VstrGDALIUCProjection[nLayer] = strProjection;
         m_VstrGDALIUCDataType[nLayer] = strDataType;
         break;
 
      case (FINE_CONS_RASTER):
         // Initial Consolidated Fine Sediment GIS data
         m_VstrGDALICFDriverCode[nLayer] = strDriverCode;
         m_VstrGDALICFDriverDesc[nLayer] = strDriverDesc;
         m_VstrGDALICFProjection[nLayer] = strProjection;
         m_VstrGDALICFDataType[nLayer] = strDataType;
         break;
 
      case (SAND_CONS_RASTER):
         // Initial Consolidated Sand Sediment GIS data
         m_VstrGDALICSDriverCode[nLayer] = strDriverCode;
         m_VstrGDALICSDriverDesc[nLayer] = strDriverDesc;
         m_VstrGDALICSProjection[nLayer] = strProjection;
         m_VstrGDALICSDataType[nLayer] = strDataType;
         break;
 
      case (COARSE_CONS_RASTER):
         // Initial Consolidated Coarse Sediment GIS data
         m_VstrGDALICCDriverCode[nLayer] = strDriverCode;
         m_VstrGDALICCDriverDesc[nLayer] = strDriverDesc;
         m_VstrGDALICCProjection[nLayer] = strProjection;
         m_VstrGDALICCDataType[nLayer] = strDataType;
         break;
      }
 
      // If present, get the missing value setting
      string const strTmp = strToLower(&strDataType);
 
      if (strTmp.find("int") != string::npos)
      {
         // This is an integer layer
         CPLPushErrorHandler(CPLQuietErrorHandler);                          // Needed to get next line to fail silently, if it fails
         m_nGISMissingValue = static_cast<int>(pGDALBand->GetNoDataValue()); // Note will fail for some formats
         CPLPopErrorHandler();
 
         if (m_nGISMissingValue != m_nMissingValue)
         {
            cerr << "   " << NOTE << "NODATA value in " << strGISFile << " is " << m_nGISMissingValue << "\n         instead using CoatalME's default integer NODATA value " << m_nMissingValue << endl;
         }
      }
 
      else
      {
         // This is a floating point layer
         CPLPushErrorHandler(CPLQuietErrorHandler);        // Needed to get next line to fail silently, if it fails
         m_dGISMissingValue = pGDALBand->GetNoDataValue(); // Note will fail for some formats
         CPLPopErrorHandler();
 
         if (! bFPIsEqual(m_dGISMissingValue, m_dMissingValue, TOLERANCE))
         {
            cerr << "   " << NOTE << "NODATA value in " << strGISFile << " is " << m_dGISMissingValue << "\n         instead using CoastalME's default floating-point NODATA value " << m_dMissingValue << endl;
         }
      }
 
      // Allocate memory for a 1D array, to hold the scan line for GDAL
      double *pdScanline = new double[m_nXGridSize];
 
      if (NULL == pdScanline)
      {
         // Error, can't allocate memory
         cerr << ERR << "cannot allocate memory for " << m_nXGridSize << " x 1D array" << endl;
         return (RTN_ERR_MEMALLOC);
      }
 
      // Now read in the data
      int nMissing = 0;
 
      for (int nY = 0; nY < m_nYGridSize; nY++)
      {
         // Read scanline
         if (CE_Failure == pGDALBand->RasterIO(GF_Read, 0, nY, m_nXGridSize, 1, pdScanline, m_nXGridSize, 1, GDT_Float64, 0, 0, NULL))
         {
            // Error while reading scanline
            cerr << ERR << CPLGetLastErrorMsg() << " in " << strGISFile << endl;
            return (RTN_ERR_RASTER_FILE_READ);
         }
 
         // All OK, so read scanline into cells (including any missing values)
         for (int nX = 0; nX < m_nXGridSize; nX++)
         {
            int nTmp;
 
            switch (nDataItem)
            {
            case (LANDFORM_RASTER):
               // Initial Landform Class GIS data, is integer TODO 030 Do we also need a landform sub-category input?
               nTmp = static_cast<int>(pdScanline[nX]);
 
               if ((isnan(nTmp)) || (nTmp == m_nGISMissingValue))
               {
                  nTmp = m_nMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->SetLFCategory(nTmp);
               break;
 
            case (INTERVENTION_CLASS_RASTER):
               // Intervention class, is integer
               nTmp = static_cast<int>(pdScanline[nX]);
 
               if ((isnan(nTmp)) || (nTmp == m_nGISMissingValue))
               {
                  nTmp = m_nMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].SetInterventionClass(nTmp);
               break;
 
            case (INTERVENTION_HEIGHT_RASTER):
               // Intervention height
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].SetInterventionHeight(dTmp);
               break;
 
            case (SUSP_SED_RASTER):
               // Initial Suspended Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].SetSuspendedSediment(dTmp);
               break;
 
            case (FINE_UNCONS_RASTER):
               // Initial Unconsolidated Fine Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetUnconsolidatedSediment()->SetFineDepth(dTmp);
               break;
 
            case (SAND_UNCONS_RASTER):
               // Initial Unconsolidated Sand Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dTmp);
               break;
 
            case (COARSE_UNCONS_RASTER):
               // Initial Unconsolidated Coarse Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dTmp);
               break;
 
            case (FINE_CONS_RASTER):
               // Initial Consolidated Fine Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetConsolidatedSediment()->SetFineDepth(dTmp);
               break;
 
            case (SAND_CONS_RASTER):
               // Initial Consolidated Sand Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetConsolidatedSediment()->SetSandDepth(dTmp);
               break;
 
            case (COARSE_CONS_RASTER):
               // Initial Consolidated Coarse Sediment GIS data
               dTmp = pdScanline[nX];
 
               if ((isnan(dTmp)) || (bFPIsEqual(dTmp, m_dGISMissingValue, TOLERANCE)))
               {
                  dTmp = m_dMissingValue;
                  nMissing++;
               }
 
               m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetConsolidatedSediment()->SetCoarseDepth(dTmp);
               break;
            }
         }
      }
 
      // Finished, so get rid of dataset object
      GDALClose(pGDALDataset);
 
      // Get rid of memory allocated to this array
      delete[] pdScanline;
 
      if (nMissing > 0)
      {
         cerr << WARN << nMissing << " missing values in " << strGISFile << endl;
         LogStream << WARN << nMissing << " missing values in " << strGISFile << endl;
      }
   }
 
   return RTN_OK;
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bWriteRasterGISFile(int const nDataItem, string const *strPlotTitle, int const nLayer, double const dElev)
{
   bool bIsInteger = false;
 
   // Begin constructing the file name for this save
   string
       strFilePathName(m_strOutPath),
       strLayer = "_layer_";
 
   stringstream ststrTmp;
 
   strLayer.append(to_string(nLayer + 1));
 
   switch (nDataItem)
   {
   case (RASTER_PLOT_BASEMENT_ELEVATION):
      strFilePathName.append(RASTER_BASEMENT_ELEVATION_NAME);
      break;
 
   case (RASTER_PLOT_SEDIMENT_TOP_ELEVATION_ELEV):
      strFilePathName.append(RASTER_SEDIMENT_TOP_NAME);
      break;
 
   case (RASTER_PLOT_OVERALL_TOP_ELEVATION):
      strFilePathName.append(RASTER_TOP_NAME);
      break;
 
   case (RASTER_PLOT_LOCAL_SLOPE_OF_CONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_LOCAL_SLOPE_NAME);
      break;
 
   case (RASTER_PLOT_SLOPE):
      strFilePathName.append(RASTER_SLOPE_NAME);
      break;
 
   case (RASTER_PLOT_CLIFF):
      strFilePathName.append(RASTER_CLIFF_NAME);
      break;
 
   case (RASTER_PLOT_SEA_DEPTH):
      strFilePathName.append(RASTER_SEA_DEPTH_NAME);
      break;
 
   case (RASTER_PLOT_AVG_SEA_DEPTH):
      strFilePathName.append(RASTER_AVG_SEA_DEPTH_NAME);
      break;
 
   case (RASTER_PLOT_WAVE_HEIGHT):
      strFilePathName.append(RASTER_WAVE_HEIGHT_NAME);
      break;
 
   case (RASTER_PLOT_AVG_WAVE_HEIGHT):
      strFilePathName.append(RASTER_AVG_WAVE_HEIGHT_NAME);
      break;
 
   case (RASTER_PLOT_WAVE_ORIENTATION):
      strFilePathName.append(RASTER_WAVE_ORIENTATION_NAME);
      break;
 
   case (RASTER_PLOT_AVG_WAVE_ORIENTATION):
      strFilePathName.append(RASTER_AVG_WAVE_ORIENTATION_NAME);
      break;
 
   case (RASTER_PLOT_BEACH_PROTECTION):
      strFilePathName.append(RASTER_BEACH_PROTECTION_NAME);
      break;
 
   case (RASTER_PLOT_POTENTIAL_PLATFORM_EROSION):
      strFilePathName.append(RASTER_POTENTIAL_PLATFORM_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_ACTUAL_PLATFORM_EROSION):
      strFilePathName.append(RASTER_ACTUAL_PLATFORM_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_POTENTIAL_PLATFORM_EROSION):
      strFilePathName.append(RASTER_TOTAL_POTENTIAL_PLATFORM_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_ACTUAL_PLATFORM_EROSION):
      strFilePathName.append(RASTER_TOTAL_ACTUAL_PLATFORM_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_POTENTIAL_BEACH_EROSION):
      strFilePathName.append(RASTER_POTENTIAL_BEACH_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_ACTUAL_BEACH_EROSION):
      strFilePathName.append(RASTER_ACTUAL_BEACH_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_POTENTIAL_BEACH_EROSION):
      strFilePathName.append(RASTER_TOTAL_POTENTIAL_BEACH_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_ACTUAL_BEACH_EROSION):
      strFilePathName.append(RASTER_TOTAL_ACTUAL_BEACH_EROSION_NAME);
      break;
 
   case (RASTER_PLOT_BEACH_DEPOSITION):
      strFilePathName.append(RASTER_BEACH_DEPOSITION_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_BEACH_DEPOSITION):
      strFilePathName.append(RASTER_TOTAL_BEACH_DEPOSITION_NAME);
      break;
 
   case (RASTER_PLOT_SUSPENDED_SEDIMENT):
      strFilePathName.append(RASTER_SUSP_SED_NAME);
      break;
 
   case (RASTER_PLOT_AVG_SUSPENDED_SEDIMENT):
      strFilePathName.append(RASTER_AVG_SUSP_SED_NAME);
      break;
 
   case (RASTER_PLOT_FINE_UNCONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_FINE_UNCONS_NAME);
      strFilePathName.append(strLayer);
      break;
 
   case (RASTER_PLOT_SAND_UNCONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_SAND_UNCONS_NAME);
      strFilePathName.append(strLayer);
      break;
 
   case (RASTER_PLOT_COARSE_UNCONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_COARSE_UNCONS_NAME);
      strFilePathName.append(strLayer);
      break;
 
   case (RASTER_PLOT_FINE_CONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_FINE_CONS_NAME);
      strFilePathName.append(strLayer);
      break;
 
   case (RASTER_PLOT_SAND_CONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_SAND_CONS_NAME);
      strFilePathName.append(strLayer);
      break;
 
   case (RASTER_PLOT_COARSE_CONSOLIDATED_SEDIMENT):
      strFilePathName.append(RASTER_COARSE_CONS_NAME);
      strFilePathName.append(strLayer);
      break;
 
   case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_FINE):
      strFilePathName.append(RASTER_CLIFF_COLLAPSE_EROSION_FINE_NAME);
      break;
 
   case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_SAND):
      strFilePathName.append(RASTER_CLIFF_COLLAPSE_EROSION_SAND_NAME);
      break;
 
   case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_COARSE):
      strFilePathName.append(RASTER_CLIFF_COLLAPSE_EROSION_COARSE_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_FINE):
      strFilePathName.append(RASTER_TOTAL_CLIFF_COLLAPSE_EROSION_FINE_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_SAND):
      strFilePathName.append(RASTER_TOTAL_CLIFF_COLLAPSE_EROSION_SAND_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_COARSE):
      strFilePathName.append(RASTER_TOTAL_CLIFF_COLLAPSE_EROSION_COARSE_NAME);
      break;
 
   case (RASTER_PLOT_CLIFF_COLLAPSE_DEPOSITION_SAND):
      strFilePathName.append(RASTER_CLIFF_COLLAPSE_DEPOSITION_SAND_NAME);
      break;
 
   case (RASTER_PLOT_CLIFF_COLLAPSE_DEPOSITION_COARSE):
      strFilePathName.append(RASTER_CLIFF_COLLAPSE_DEPOSITION_COARSE_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_DEPOSITION_SAND):
      strFilePathName.append(RASTER_TOTAL_CLIFF_COLLAPSE_DEPOSITION_SAND_NAME);
      break;
 
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_DEPOSITION_COARSE):
      strFilePathName.append(RASTER_TOTAL_CLIFF_COLLAPSE_DEPOSITION_COARSE_NAME);
      break;
 
   case (RASTER_PLOT_INTERVENTION_HEIGHT):
      strFilePathName.append(RASTER_INTERVENTION_HEIGHT_NAME);
      break;
 
   case (RASTER_PLOT_DEEP_WATER_WAVE_ORIENTATION):
      strFilePathName.append(RASTER_DEEP_WATER_WAVE_ORIENTATION_NAME);
      break;
 
   case (RASTER_PLOT_DEEP_WATER_WAVE_HEIGHT):
      strFilePathName.append(RASTER_DEEP_WATER_WAVE_HEIGHT_NAME);
      break;
 
   case (RASTER_PLOT_POLYGON_GAIN_OR_LOSS):
      strFilePathName.append(RASTER_POLYGON_GAIN_OR_LOSS_NAME);
      break;
 
   case (RASTER_PLOT_DEEP_WATER_WAVE_PERIOD):
      strFilePathName.append(RASTER_WAVE_PERIOD_NAME);
      break;
 
   case (RASTER_PLOT_SEDIMENT_INPUT):
      strFilePathName.append(RASTER_SEDIMENT_INPUT_EVENT_NAME);
      break;
 
   case (RASTER_PLOT_BEACH_MASK):
      bIsInteger = true;
      strFilePathName.append(RASTER_BEACH_MASK_NAME);
      break;
 
   case (RASTER_PLOT_POTENTIAL_PLATFORM_EROSION_MASK):
      bIsInteger = true;
      strFilePathName.append(RASTER_POTENTIAL_PLATFORM_EROSION_MASK_NAME);
      break;
 
   case (RASTER_PLOT_INUNDATION_MASK):
      bIsInteger = true;
      strFilePathName.append(RASTER_INUNDATION_MASK_NAME);
      break;
 
   case (RASTER_PLOT_SLICE):
      bIsInteger = true;
      ststrTmp.str("");
      ststrTmp.clear();
 
      // TODO 031 Get working for multiple slices
      strFilePathName.append(RASTER_SLICE_NAME);
      ststrTmp << "_" << dElev << "_";
      strFilePathName.append(ststrTmp.str());
      break;
 
   case (RASTER_PLOT_LANDFORM):
      bIsInteger = true;
      strFilePathName.append(RASTER_LANDFORM_NAME);
      break;
 
   case (RASTER_PLOT_INTERVENTION_CLASS):
      bIsInteger = true;
      strFilePathName.append(RASTER_INTERVENTION_CLASS_NAME);
      break;
 
   case (RASTER_PLOT_COAST):
      bIsInteger = true;
      strFilePathName.append(RASTER_COAST_NAME);
      break;
 
   case (RASTER_PLOT_NORMAL_PROFILE):
      bIsInteger = true;
      strFilePathName.append(RASTER_COAST_NORMAL_NAME);
      break;
 
   case (RASTER_PLOT_ACTIVE_ZONE):
      bIsInteger = true;
      strFilePathName.append(RASTER_ACTIVE_ZONE_NAME);
      break;
 
   case (RASTER_PLOT_POLYGON):
      bIsInteger = true;
      strFilePathName.append(RASTER_POLYGON_NAME);
      break;
 
   case (RASTER_PLOT_SHADOW_ZONE):
      bIsInteger = true;
      strFilePathName.append(RASTER_SHADOW_ZONE_NAME);
      break;
 
   case (RASTER_PLOT_SHADOW_DOWNDRIFT_ZONE):
      bIsInteger = true;
      strFilePathName.append(RASTER_SHADOW_DOWNDRIFT_ZONE_NAME);
      break;
 
   case (RASTER_PLOT_POLYGON_UPDRIFT_OR_DOWNDRIFT):
      bIsInteger = true;
      strFilePathName.append(RASTER_POLYGON_UPDRIFT_OR_DOWNDRIFT_NAME);
      break;
 
   case (RASTER_PLOT_SETUP_SURGE_FLOOD_MASK):
      bIsInteger = true;
      strFilePathName.append(RASTER_SETUP_SURGE_FLOOD_MASK_NAME);
      break;
 
   case (RASTER_PLOT_SETUP_SURGE_RUNUP_FLOOD_MASK):
      bIsInteger = true;
      strFilePathName.append(RASTER_SETUP_SURGE_RUNUP_FLOOD_MASK_NAME);
      break;
 
   case (RASTER_PLOT_WAVE_FLOOD_LINE):
      bIsInteger = true;
      strFilePathName.append(RASTER_WAVE_FLOOD_LINE_NAME);
      break;
   }
 
   // Append the 'save number' to the filename, and prepend zeros to the save number
   ststrTmp.str("");
   ststrTmp.clear();
 
   strFilePathName.append("_");
 
   if (m_bGISSaveDigitsSequential)
   {
      // Save number is m_bGISSaveDigitsSequential
      ststrTmp << FillToWidth('0', m_nGISMaxSaveDigits) << m_nGISSave;
   }
 
   else
   {
      // Save number is iteration
      ststrTmp << FillToWidth('0', m_nGISMaxSaveDigits) << m_ulIter;
   }
 
   strFilePathName.append(ststrTmp.str());
 
   // Finally, maybe append the extension
   if (! m_strGDALRasterOutputDriverExtension.empty())
   {
      strFilePathName.append(".");
      strFilePathName.append(m_strGDALRasterOutputDriverExtension);
   }
 
   // TODO 065 Used to try to debug floating point exception in pDriver->Create() below
   // CPLSetConfigOption("CPL_DEBUG", "ON");
   // CPLSetConfigOption("GDAL_NUM_THREADS", "1");
 
   GDALDriver *pDriver;
   GDALDataset *pDataSet;
 
   if (m_bGDALCanCreate)
   {
      // The user-requested raster driver supports the Create() method
      pDriver = GetGDALDriverManager()->GetDriverByName(m_strRasterGISOutFormat.c_str());
 
      if ((nDataItem == RASTER_PLOT_INUNDATION_MASK) || (nDataItem == RASTER_PLOT_SETUP_SURGE_FLOOD_MASK) || (nDataItem == RASTER_PLOT_SETUP_SURGE_RUNUP_FLOOD_MASK))
      {
         pDataSet = pDriver->Create(strFilePathName.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Int16, m_papszGDALRasterOptions);
      }
 
      else if (m_strRasterGISOutFormat == "gpkg")
      {
         // TODO 065 Floating point exception here
         pDataSet = pDriver->Create(strFilePathName.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Byte, m_papszGDALRasterOptions);
      }
 
      else
      {
         pDataSet = pDriver->Create(strFilePathName.c_str(), m_nXGridSize, m_nYGridSize, 1, m_GDALWriteFloatDataType, m_papszGDALRasterOptions);
      }
 
      if (NULL == pDataSet)
      {
         // Error, couldn't create file
         cerr << ERR << "cannot create " << m_strRasterGISOutFormat << " file named " << strFilePathName << endl
              << CPLGetLastErrorMsg() << endl;
         return false;
      }
   }
 
   else
   {
      // The user-requested raster driver does not support the Create() method, so we must first create a memory-file dataset
      pDriver = GetGDALDriverManager()->GetDriverByName("MEM");
      pDataSet = pDriver->Create("", m_nXGridSize, m_nYGridSize, 1, m_GDALWriteFloatDataType, NULL);
 
      if (NULL == pDataSet)
      {
         // Couldn't create in-memory file dataset
         cerr << ERR << "cannot create in-memory file for " << m_strRasterGISOutFormat << " file named " << strFilePathName << endl
              << CPLGetLastErrorMsg() << endl;
         return false;
      }
   }
 
   // Set projection info for output dataset (will be same as was read in from basement DEM)
   CPLPushErrorHandler(CPLQuietErrorHandler);                       // Needed to get next line to fail silently, if it fails
   pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str()); // Will fail for some formats
   CPLPopErrorHandler();
 
   // Set geotransformation info for output dataset (will be same as was read in from DEM)
   if (CE_Failure == pDataSet->SetGeoTransform(m_dGeoTransform))
      LogStream << WARN << "cannot write geotransformation information to " << m_strRasterGISOutFormat << " file named " << strFilePathName << endl
                << CPLGetLastErrorMsg() << endl;
 
   // Allocate memory for a 1D array, to hold the floating point raster band data for GDAL
   double *pdRaster = new double[m_ulNumCells];
 
   if (NULL == pdRaster)
   {
      // Error, can't allocate memory
      cerr << ERR << "cannot allocate memory for " << m_ulNumCells << " x 1D floating-point array for " << m_strRasterGISOutFormat << " file named " << strFilePathName << endl;
      return (RTN_ERR_MEMALLOC);
   }
 
   bool bScaleOutput = false;
   double dRangeScale = 0;
   double dDataMin = 0;
 
   if (! m_bGDALCanWriteFloat)
   {
      double dDataMax = 0;
 
      // The output file format cannot handle floating-point numbers, so we may need to scale the output
      GetRasterOutputMinMax(nDataItem, dDataMin, dDataMax, nLayer, 0);
 
      double const dDataRange = dDataMax - dDataMin;
      double const dWriteRange = static_cast<double>(m_lGDALMaxCanWrite - m_lGDALMinCanWrite);
 
      if (dDataRange > 0)
         dRangeScale = dWriteRange / dDataRange;
 
      // If we are attempting to write values which are outside this format's allowable range, and the user has set the option, then scale the output
      if (((dDataMin < static_cast<double>(m_lGDALMinCanWrite)) || (dDataMax > static_cast<double>(m_lGDALMaxCanWrite))) && m_bScaleRasterOutput)
         bScaleOutput = true;
   }
 
   // Fill the array
   int
       n = 0,
       nPoly = 0,
       nTopLayer = 0;
   double dTmp = 0;
 
   for (int nY = 0; nY < m_nYGridSize; nY++)
   {
      for (int nX = 0; nX < m_nXGridSize; nX++)
      {
         switch (nDataItem)
         {
         case (RASTER_PLOT_BASEMENT_ELEVATION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetBasementElev();
            break;
 
         case (RASTER_PLOT_SEDIMENT_TOP_ELEVATION_ELEV):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
            break;
 
         case (RASTER_PLOT_OVERALL_TOP_ELEVATION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetSedimentPlusInterventionTopElev();
            break;
 
         case (RASTER_PLOT_LOCAL_SLOPE_OF_CONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetLocalConsSlope();
            break;
 
         case (RASTER_PLOT_SLOPE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetSlope();
            break;
 
         case (RASTER_PLOT_CLIFF):
            dTmp = static_cast<double>(m_pRasterGrid->m_Cell[nX][nY].bIsCliff());
            break;
 
         case (RASTER_PLOT_SEA_DEPTH):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetSeaDepth();
            break;
 
         case (RASTER_PLOT_AVG_SEA_DEPTH):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotSeaDepth() / static_cast<double>(m_ulIter);
            break;
 
         case (RASTER_PLOT_WAVE_HEIGHT):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_AVG_WAVE_HEIGHT):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotWaveHeight() / static_cast<double>(m_ulIter);
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_WAVE_ORIENTATION):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_AVG_WAVE_ORIENTATION):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotWaveAngle() / static_cast<double>(m_ulIter);
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_BEACH_PROTECTION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetBeachProtectionFactor();
 
            if (bFPIsEqual(dTmp, DBL_NODATA, TOLERANCE))
               dTmp = m_dMissingValue;
 
            else
               dTmp = 1 - dTmp; // Output the inverse, seems more intuitive
 
            break;
 
         case (RASTER_PLOT_POTENTIAL_PLATFORM_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetPotentialPlatformErosion();
            break;
 
         case (RASTER_PLOT_ACTUAL_PLATFORM_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetActualPlatformErosion();
            break;
 
         case (RASTER_PLOT_TOTAL_POTENTIAL_PLATFORM_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotPotentialPlatformErosion();
            break;
 
         case (RASTER_PLOT_TOTAL_ACTUAL_PLATFORM_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotActualPlatformErosion();
            break;
 
         case (RASTER_PLOT_POTENTIAL_BEACH_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetPotentialBeachErosion();
            break;
 
         case (RASTER_PLOT_ACTUAL_BEACH_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetActualBeachErosion();
            break;
 
         case (RASTER_PLOT_TOTAL_POTENTIAL_BEACH_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotPotentialBeachErosion();
            break;
 
         case (RASTER_PLOT_TOTAL_ACTUAL_BEACH_EROSION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotActualBeachErosion();
            break;
 
         case (RASTER_PLOT_BEACH_DEPOSITION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetBeachDeposition();
            break;
 
         case (RASTER_PLOT_TOTAL_BEACH_DEPOSITION):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotBeachDeposition();
            break;
 
         case (RASTER_PLOT_SUSPENDED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetSuspendedSediment();
            break;
 
         case (RASTER_PLOT_AVG_SUSPENDED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotSuspendedSediment() / static_cast<double>(m_ulIter);
            break;
 
         case (RASTER_PLOT_FINE_UNCONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetUnconsolidatedSediment()->dGetFineDepth();
            break;
 
         case (RASTER_PLOT_SAND_UNCONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
            break;
 
         case (RASTER_PLOT_COARSE_UNCONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
            break;
 
         case (RASTER_PLOT_FINE_CONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetConsolidatedSediment()->dGetFineDepth();
            break;
 
         case (RASTER_PLOT_SAND_CONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetConsolidatedSediment()->dGetSandDepth();
            break;
 
         case (RASTER_PLOT_COARSE_CONSOLIDATED_SEDIMENT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nLayer)->pGetConsolidatedSediment()->dGetCoarseDepth();
            break;
 
         case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_FINE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseErosionFine();
            break;
 
         case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_SAND):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseErosionSand();
            break;
 
         case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_COARSE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseErosionCoarse();
            break;
 
         case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_FINE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotCliffCollapseFine();
            break;
 
         case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_SAND):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotCliffCollapseSand();
            break;
 
         case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_COARSE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotCliffCollapseCoarse();
            break;
 
         case (RASTER_PLOT_CLIFF_COLLAPSE_DEPOSITION_SAND):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseSandTalusDeposition();
            break;
 
         case (RASTER_PLOT_CLIFF_COLLAPSE_DEPOSITION_COARSE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetThisIterCliffCollapseCoarseTalusDeposition();
            break;
 
         case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_DEPOSITION_SAND):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotSandTalusDeposition();
            break;
 
         case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_DEPOSITION_COARSE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetTotCoarseTalusDeposition();
            break;
 
         case (RASTER_PLOT_INTERVENTION_HEIGHT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetInterventionHeight();
            break;
 
         case (RASTER_PLOT_DEEP_WATER_WAVE_ORIENTATION):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveAngle();
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_DEEP_WATER_WAVE_HEIGHT):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveHeight();
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_DEEP_WATER_WAVE_PERIOD):
            if (m_pRasterGrid->m_Cell[nX][nY].bIsInundated())
               dTmp = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWavePeriod();
 
            else
               dTmp = 0;
 
            break;
 
         case (RASTER_PLOT_POLYGON_GAIN_OR_LOSS):
            nPoly = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
 
            if (nPoly == INT_NODATA)
               dTmp = m_dMissingValue;
 
            else
            {
               // Get total volume (all sediment size classes) of change in sediment for this polygon for this timestep (-ve erosion, +ve deposition)
               dTmp = m_pVCoastPolygon[nPoly]->dGetBeachDepositionAndSuspensionAllUncons() * m_dCellArea;
 
               // Calculate the rate in m^3 / sec
               dTmp /= (m_dTimeStep * 3600);
            }
 
            break;
 
         case (RASTER_PLOT_POTENTIAL_PLATFORM_EROSION_MASK):
            // cppcheck-suppress assignBoolToFloat
            dTmp = m_pRasterGrid->m_Cell[nX][nY].bPotentialPlatformErosion();
            break;
 
         case (RASTER_PLOT_INUNDATION_MASK):
            // cppcheck-suppress assignBoolToFloat
            dTmp = m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea();
            break;
 
         case (RASTER_PLOT_BEACH_MASK):
            dTmp = 0;
            nTopLayer = m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
 
            if ((nTopLayer == INT_NODATA) || (nTopLayer == NO_NONZERO_THICKNESS_LAYERS))
               break;
 
            if ((m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->dGetUnconsolidatedThickness() > 0) && (m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev() > m_dThisIterSWL))
               dTmp = 1;
 
            break;
 
         case (RASTER_PLOT_SLICE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].nGetLayerAtElev(dElev);
            break;
 
         case (RASTER_PLOT_LANDFORM):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->nGetLFCategory();
 
            if ((static_cast<int>(dTmp) == LF_CAT_DRIFT) || (static_cast<int>(dTmp) == LF_CAT_CLIFF))
               dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->nGetLFSubCategory();
 
            break;
 
         case (RASTER_PLOT_INTERVENTION_CLASS):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].nGetInterventionClass();
            break;
 
         case (RASTER_PLOT_COAST):
            dTmp = (m_pRasterGrid->m_Cell[nX][nY].bIsCoastline() ? 1 : 0);
            break;
 
         case (RASTER_PLOT_NORMAL_PROFILE):
            // dTmp = (m_pRasterGrid->m_Cell[nX][nY].bIsProfile() ? 1 : 0);
            dTmp = m_pRasterGrid->m_Cell[nX][nY].nGetProfileID();
            break;
 
         case (RASTER_PLOT_ACTIVE_ZONE):
            dTmp = (m_pRasterGrid->m_Cell[nX][nY].bIsInActiveZone() ? 1 : 0);
            break;
 
         case (RASTER_PLOT_POLYGON):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
            break;
 
         case (RASTER_PLOT_SHADOW_ZONE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].nGetShadowZoneNumber();
            break;
 
         case (RASTER_PLOT_SHADOW_DOWNDRIFT_ZONE):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].nGetDownDriftZoneNumber();
            break;
 
         case (RASTER_PLOT_POLYGON_UPDRIFT_OR_DOWNDRIFT):
            nPoly = m_pRasterGrid->m_Cell[nX][nY].nGetPolygonID();
 
            if (nPoly == INT_NODATA)
               dTmp = m_nMissingValue;
 
            else
            {
               if (m_pVCoastPolygon[nPoly]->bDownCoastThisIter())
                  dTmp = 1;
 
               else
                  dTmp = 0;
            }
 
            break;
 
         case (RASTER_PLOT_SEDIMENT_INPUT):
            dTmp = m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetTotAllSedimentInputDepth();
            break;
 
         case (RASTER_PLOT_SETUP_SURGE_FLOOD_MASK):
            dTmp = (m_pRasterGrid->m_Cell[nX][nY].bIsFloodBySetupSurge() ? 1 : 0);
            break;
 
         case (RASTER_PLOT_SETUP_SURGE_RUNUP_FLOOD_MASK):
            dTmp = (m_pRasterGrid->m_Cell[nX][nY].bIsFloodBySetupSurgeRunup() ? 1 : 0);
            break;
 
         case (RASTER_PLOT_WAVE_FLOOD_LINE):
            dTmp = (m_pRasterGrid->m_Cell[nX][nY].bIsFloodline() ? 1 : 0);
            break;
         }
 
         // If necessary, scale this value
         if (bScaleOutput)
         {
            if (bFPIsEqual(dTmp, DBL_NODATA, TOLERANCE))
               dTmp = 0; // TODO 032 Improve this
 
            else
               dTmp = dRound(static_cast<double>(m_lGDALMinCanWrite) + (dRangeScale * (dTmp - dDataMin)));
         }
 
         // Write this value to the array
         pdRaster[n++] = dTmp;
      }
   }
 
   // Create a single raster band
   GDALRasterBand *pBand = pDataSet->GetRasterBand(1);
 
   // And fill it with the NODATA value
   if (bIsInteger)
      pBand->Fill(m_nMissingValue);
 
   else
      pBand->Fill(m_dMissingValue);
 
   // Set value units for this band
   string strUnits;
 
   switch (nDataItem)
   {
   case (RASTER_PLOT_BASEMENT_ELEVATION):
   case (RASTER_PLOT_SEDIMENT_TOP_ELEVATION_ELEV):
   case (RASTER_PLOT_OVERALL_TOP_ELEVATION):
   case (RASTER_PLOT_SEA_DEPTH):
   case (RASTER_PLOT_AVG_SEA_DEPTH):
   case (RASTER_PLOT_WAVE_HEIGHT):
   case (RASTER_PLOT_AVG_WAVE_HEIGHT):
   case (RASTER_PLOT_POTENTIAL_PLATFORM_EROSION):
   case (RASTER_PLOT_ACTUAL_PLATFORM_EROSION):
   case (RASTER_PLOT_TOTAL_POTENTIAL_PLATFORM_EROSION):
   case (RASTER_PLOT_TOTAL_ACTUAL_PLATFORM_EROSION):
   case (RASTER_PLOT_POTENTIAL_BEACH_EROSION):
   case (RASTER_PLOT_ACTUAL_BEACH_EROSION):
   case (RASTER_PLOT_TOTAL_POTENTIAL_BEACH_EROSION):
   case (RASTER_PLOT_TOTAL_ACTUAL_BEACH_EROSION):
   case (RASTER_PLOT_BEACH_DEPOSITION):
   case (RASTER_PLOT_TOTAL_BEACH_DEPOSITION):
   case (RASTER_PLOT_SUSPENDED_SEDIMENT):
   case (RASTER_PLOT_AVG_SUSPENDED_SEDIMENT):
   case (RASTER_PLOT_FINE_UNCONSOLIDATED_SEDIMENT):
   case (RASTER_PLOT_SAND_UNCONSOLIDATED_SEDIMENT):
   case (RASTER_PLOT_COARSE_UNCONSOLIDATED_SEDIMENT):
   case (RASTER_PLOT_FINE_CONSOLIDATED_SEDIMENT):
   case (RASTER_PLOT_SAND_CONSOLIDATED_SEDIMENT):
   case (RASTER_PLOT_COARSE_CONSOLIDATED_SEDIMENT):
   case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_FINE):
   case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_SAND):
   case (RASTER_PLOT_CLIFF_COLLAPSE_EROSION_COARSE):
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_FINE):
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_SAND):
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_EROSION_COARSE):
   case (RASTER_PLOT_CLIFF_COLLAPSE_DEPOSITION_SAND):
   case (RASTER_PLOT_CLIFF_COLLAPSE_DEPOSITION_COARSE):
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_DEPOSITION_SAND):
   case (RASTER_PLOT_TOTAL_CLIFF_COLLAPSE_DEPOSITION_COARSE):
   case (RASTER_PLOT_INTERVENTION_HEIGHT):
   case (RASTER_PLOT_DEEP_WATER_WAVE_HEIGHT):
   case (RASTER_PLOT_SEDIMENT_INPUT):
      strUnits = "m";
      break;
 
   case (RASTER_PLOT_LOCAL_SLOPE_OF_CONSOLIDATED_SEDIMENT):
      strUnits = "m/m";
      break;
 
   case (RASTER_PLOT_WAVE_ORIENTATION):
   case (RASTER_PLOT_AVG_WAVE_ORIENTATION):
      strUnits = "degrees";
      break;
 
   case (RASTER_PLOT_POLYGON_GAIN_OR_LOSS):
      strUnits = "cumecs";
      break;
 
   case (RASTER_PLOT_DEEP_WATER_WAVE_PERIOD):
      strUnits = "secs";
      break;
 
   case (RASTER_PLOT_POTENTIAL_PLATFORM_EROSION_MASK):
   case (RASTER_PLOT_INUNDATION_MASK):
   case (RASTER_PLOT_BEACH_MASK):
   case (RASTER_PLOT_SLICE):
   case (RASTER_PLOT_LANDFORM):
   case (RASTER_PLOT_INTERVENTION_CLASS):
   case (RASTER_PLOT_COAST):
   case (RASTER_PLOT_NORMAL_PROFILE):
   case (RASTER_PLOT_ACTIVE_ZONE):
   case (RASTER_PLOT_POLYGON):
   case (RASTER_PLOT_SHADOW_ZONE):
   case (RASTER_PLOT_SHADOW_DOWNDRIFT_ZONE):
   case (RASTER_PLOT_POLYGON_UPDRIFT_OR_DOWNDRIFT):
   case (RASTER_PLOT_SETUP_SURGE_FLOOD_MASK):
   case (RASTER_PLOT_SETUP_SURGE_RUNUP_FLOOD_MASK):
   case (RASTER_PLOT_WAVE_FLOOD_LINE):
      strUnits = "none";
      break;
   }
 
   CPLPushErrorHandler(CPLQuietErrorHandler); // Needed to get next line to fail silently, if it fails
   pBand->SetUnitType(strUnits.c_str());      // Not supported for some GIS formats
   CPLPopErrorHandler();
 
   // Tell the output dataset about NODATA (missing values)
   CPLPushErrorHandler(CPLQuietErrorHandler); // Needed to get next line to fail silently, if it fails
 
   if (bIsInteger)
      pBand->SetNoDataValue(m_nMissingValue); // Will fail for some formats
 
   else
      pBand->SetNoDataValue(m_dMissingValue); // Will fail for some formats
 
   CPLPopErrorHandler();
 
   // Construct the description
   string strDesc(*strPlotTitle);
 
   if (nDataItem == RASTER_PLOT_SLICE)
   {
      ststrTmp.clear();
      ststrTmp << dElev << "m, ";
      strDesc.append(ststrTmp.str());
   }
 
   strDesc.append(" at ");
   strDesc.append(strDispTime(m_dSimElapsed, false, false));
 
   // Set the GDAL description
   pBand->SetDescription(strDesc.c_str());
 
   // Set raster category names
   char **papszCategoryNames = NULL;
 
   switch (nDataItem)
   {
   case (RASTER_PLOT_SLICE):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Basement");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 0");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 1");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 2");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 3");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 4");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 5");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 6");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 7");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 8");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Layer 9");
      break;
 
   case (RASTER_PLOT_LANDFORM):
      papszCategoryNames = CSLAddString(papszCategoryNames, "None");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Hinterland");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Sea");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Cliff");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Drift");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Intervention");
 
      papszCategoryNames = CSLAddString(papszCategoryNames, "Cliff on Coastline");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Inland Cliff");
 
      papszCategoryNames = CSLAddString(papszCategoryNames, "Mixed Drift");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Talus");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Beach");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Dunes");
      break;
 
   case (RASTER_PLOT_INTERVENTION_CLASS):
      papszCategoryNames = CSLAddString(papszCategoryNames, "None");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Structural");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Non-Structural");
      break;
 
   case (RASTER_PLOT_COAST):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not coastline");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Coastline");
      break;
 
   case (RASTER_PLOT_NORMAL_PROFILE):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not coastline-normal profile");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Coastline-normal profile");
      break;
 
   case (RASTER_PLOT_ACTIVE_ZONE):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not in active zone");
      papszCategoryNames = CSLAddString(papszCategoryNames, "In active zone");
      break;
 
   case (RASTER_PLOT_POLYGON):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not polygon");
      papszCategoryNames = CSLAddString(papszCategoryNames, "In polygon");
      break;
 
   case (RASTER_PLOT_SHADOW_ZONE):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not in shadow zone");
      papszCategoryNames = CSLAddString(papszCategoryNames, "In shadow zone");
      break;
 
   case (RASTER_PLOT_SHADOW_DOWNDRIFT_ZONE):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not in shadow downdrift zone");
      papszCategoryNames = CSLAddString(papszCategoryNames, "In shadow downdrift zone");
      break;
 
   case (RASTER_PLOT_POLYGON_UPDRIFT_OR_DOWNDRIFT):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Updrift movement of unconsolidated sediment ");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Downdrift movement of unconsolidated sediment");
      break;
 
   case (RASTER_PLOT_SETUP_SURGE_FLOOD_MASK):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Inundated by swl setup and surge ");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not inundated by swl setup and surge");
      break;
 
   case (RASTER_PLOT_SETUP_SURGE_RUNUP_FLOOD_MASK):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Inundated by swl setup, surge and runup ");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not inundated by swl setup, surge and runup");
      break;
 
   case (RASTER_PLOT_WAVE_FLOOD_LINE):
      papszCategoryNames = CSLAddString(papszCategoryNames, "Intersection line of inundation ");
      papszCategoryNames = CSLAddString(papszCategoryNames, "Not inundated by swl waves and runup");
      break;
   }
 
   CPLPushErrorHandler(CPLQuietErrorHandler);   // Needed to get next line to fail silently, if it fails
   pBand->SetCategoryNames(papszCategoryNames); // Not supported for some GIS formats
   CPLPopErrorHandler();
 
   // Now write the data with optimized I/O
   // Enable multi-threaded compression for faster writing
   CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
 
   if (CE_Failure == pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL))
   {
      // Write error, better error message
      cerr << ERR << "cannot write data for " << m_strRasterGISOutFormat << " file named " << strFilePathName << endl
           << CPLGetLastErrorMsg() << endl;
      delete[] pdRaster;
      return false;
   }
 
   // Calculate statistics for this band
   double dMin, dMax, dMean, dStdDev;
   CPLPushErrorHandler(CPLQuietErrorHandler); // Needed to get next line to fail silently, if it fails
   pBand->ComputeStatistics(false, &dMin, &dMax, &dMean, &dStdDev, NULL, NULL);
   CPLPopErrorHandler();
 
   // And then write the statistics
   CPLPushErrorHandler(CPLQuietErrorHandler); // Needed to get next line to fail silently, if it fails
   pBand->SetStatistics(dMin, dMax, dMean, dStdDev);
   CPLPopErrorHandler();
 
   if (! m_bGDALCanCreate)
   {
      // Since the user-selected raster driver cannot use the Create() method, we have been writing to a dataset created by the in-memory driver. So now we need to use CreateCopy() to copy this in-memory dataset to a file in the user-specified raster driver format
      GDALDriver *pOutDriver = GetGDALDriverManager()->GetDriverByName(m_strRasterGISOutFormat.c_str());
      GDALDataset *pOutDataSet = pOutDriver->CreateCopy(strFilePathName.c_str(), pDataSet, false, m_papszGDALRasterOptions, NULL, NULL);
 
      if (NULL == pOutDataSet)
      {
         // Couldn't create file
         cerr << ERR << "cannot create " << m_strRasterGISOutFormat << " file named " << strFilePathName << endl
              << CPLGetLastErrorMsg() << endl;
         return false;
      }
 
      // Get rid of this user-selected dataset object
      GDALClose(pOutDataSet);
   }
 
   // Get rid of dataset object
   GDALClose(pDataSet);
 
   // Also get rid of memory allocated to this array
   delete[] pdRaster;
 
   return true;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nInterpolateWavesToPolygonCells(vector<double> const *pVdX, vector<double> const *pVdY, vector<double> const *pVdHeightX, vector<double> const *pVdHeightY)
{
   int nXSize = 0;
   int nYSize = 0;
 
   double dXAvg = 0;
   double dYAvg = 0;
 
   nXSize = m_nXMaxBoundingBox - m_nXMinBoundingBox + 1;
   nYSize = m_nYMaxBoundingBox - m_nYMinBoundingBox + 1;
   int const nGridSize = nXSize * nYSize;
 
   unsigned int const nPoints = static_cast<unsigned int>(pVdX->size());
 
   //    // DEBUG CODE ============================================================================================================
   // for (int nn = 0; nn < nPoints; nn++)
   // {
   // LogStream << nn << " " << dX[nn] << " " << dY[nn] << " " << dZ[nn] << endl;
   // }
   //
   // m_nXMaxBoundingBox = m_nXGridSize-1;
   // m_nYMaxBoundingBox = m_nYGridSize-1;
   // m_nXMinBoundingBox = 0;
   // m_nYMinBoundingBox = 0;
   //    // DEBUG CODE ============================================================================================================
 
   vector<double> VdOutX(nGridSize, 0);
   vector<double> VdOutY(nGridSize, 0);
 
   // Do first for X, then for Y
   for (int nDirection = 0; nDirection < 2; nDirection++)
   {
      // Use the GDALGridCreate() linear interpolation algorithm: this computes a Delaunay triangulation of the point cloud, finding in which triangle of the triangulation the point is, and by doing linear interpolation from its barycentric coordinates within the triangle. If the point is not in any triangle, depending on the radius, the algorithm will use the value of the nearest point or the nodata value. Only available in GDAL 2.1 and later TODO 086
 
      // Performance optimization: Use optimized settings for large grids
      GDALGridLinearOptions *pOptions = new GDALGridLinearOptions();
      pOptions->dfNoDataValue = m_dMissingValue;                  // Set the no-data marker to fill empty points
      pOptions->dfRadius = -1;                                    // Set the search radius to infinite
      pOptions->nSizeOfStructure = sizeof(GDALGridLinearOptions); // Needed for GDAL 3.6 onwards, see https://gdal.org/api/gdal_alg.html#_CPPv421GDALGridLinearOptions
 
      // Performance enhancement: Enable grid threading for this specific operation
      CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
 
      // pOptions.dfRadius = static_cast<double>(nXSize + nYSize) / 2.0;                       // Set the search radius
 
      // GDALGridNearestNeighborOptions* pOptions = new GDALGridNearestNeighborOptions();
      // pOptions->dfNoDataValue = m_dMissingValue; // Set the no-data marker to fill empty points
      // pOptions->dfRadius = -1;                   // Set the search radius to infinite
 
      // Call GDALGridCreate() TODO 086
      int nRet;
 
      if (nDirection == 0)
      {
         nRet = GDALGridCreate(GGA_Linear, pOptions, nPoints, pVdX->data(), pVdY->data(), pVdHeightX->data(), m_nXMinBoundingBox, m_nXMaxBoundingBox, m_nYMinBoundingBox, m_nYMaxBoundingBox, nXSize, nYSize, GDT_Float64, VdOutX.data(), NULL, NULL);
      }
 
      else
      {
         nRet = GDALGridCreate(GGA_Linear, pOptions, nPoints, pVdX->data(), pVdY->data(), pVdHeightY->data(), m_nXMinBoundingBox, m_nXMaxBoundingBox, m_nYMinBoundingBox, m_nYMaxBoundingBox, nXSize, nYSize, GDT_Float64, VdOutY.data(), NULL, NULL);
      }
 
      // int nRet;
      // if (nDirection == 0)
      // {
      // nRet = GDALGridCreate(GGA_NearestNeighbor, pOptions, nPoints, pVdX->data(), pVdY->data(), pVdHeightX->data(), m_nXMinBoundingBox, m_nXMaxBoundingBox, m_nYMinBoundingBox, m_nYMaxBoundingBox, nXSize, nYSize, GDT_Float64, VdOutX.data(), NULL, NULL);
      // }
      // else
      // {
      // nRet = GDALGridCreate(GGA_NearestNeighbor, pOptions, nPoints, pVdX->data(), pVdY->data(), pVdHeightY->data(), m_nXMinBoundingBox, m_nXMaxBoundingBox, m_nYMinBoundingBox, m_nYMaxBoundingBox, nXSize, nYSize, GDT_Float64, VdOutY.data(), NULL, NULL);
      // }
 
      delete pOptions;
 
      if (nRet == CE_Failure)
      {
         cerr << CPLGetLastErrorMsg() << endl;
         return RTN_ERR_GRIDCREATE;
      }
 
      if (nDirection == 0)
      {
         int nXValid = 0;
 
         // Safety check: unfortunately, GDALGridCreate(() outputs NaNs and other crazy values when the polygons are far from regular. So check for these
         for (unsigned int n = 0; n < VdOutX.size(); n++)
         {
            if (isnan(VdOutX[n]))
               VdOutX[n] = m_dMissingValue;
 
            else if (tAbs(VdOutX[n]) > 1e10)
               VdOutX[n] = m_dMissingValue;
 
            else
            {
               dXAvg += VdOutX[n];
               nXValid++;
            }
         }
 
         dXAvg /= nXValid;
      }
 
      else
      {
         int nYValid = 0;
 
         // Safety check: unfortunately, GDALGridCreate(() outputs NaNs and other crazy values when the polygon are far from regular. So check for these
         for (unsigned int n = 0; n < VdOutY.size(); n++)
         {
            if (isnan(VdOutY[n]))
               VdOutY[n] = m_dMissingValue;
 
            else if (tAbs(VdOutY[n]) > 1e10)
               VdOutY[n] = m_dMissingValue;
 
            else
            {
               dYAvg += VdOutY[n];
               nYValid++;
            }
         }
 
         dYAvg /= nYValid;
      }
 
      // // DEBUG CODE ===========================================================================================================
      // string strOutFile = m_strOutPath;
      // strOutFile += "sea_wave_interpolation_";
      // if (nDirection == 0)
      // strOutFile += "X_";
      // else
      // strOutFile += "Y_";
      // strOutFile += to_string(m_ulIter);
      // strOutFile += ".tif";
      //
      // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
      // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
      // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
      // pDataSet->SetGeoTransform(m_dGeoTransform);
      // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
      // int m = 0;
      // int n = 0;
      // for (int nY = 0; nY < m_nYGridSize; nY++)
      // {
      // for (int nX = 0; nX < m_nXGridSize; nX++)
      // {
      // if ((nX < m_nXMinBoundingBox) || (nY < m_nYMinBoundingBox))
      // {
      // pdRaster[n++] = DBL_NODATA;
      // }
      // else
      // {
      //          // Write this value to the array
      // if (nDirection == 0)
      // {
      // if (m < static_cast<int>(VdOutX.size()))
      // {
      // pdRaster[n++] = VdOutX[m++];
      //                // LogStream << "nDirection = " << nDirection << " [" << nX << "][" << nY << "] = " << VpdOutX[n] << endl;
      // }
      // }
      // else
      // {
      // if (m < static_cast<int>(VdOutY.size()))
      // {
      // pdRaster[n++] = VdOutY[m++];
      //                // LogStream << "nDirection = " << nDirection << " [" << nX << "][" << nY << "] = " << VpdOutY[n] << endl;
      // }
      // }
      // }
      // }
      // }
      //
      // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
      // pBand->SetNoDataValue(m_dMissingValue);
      // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
      //
      // if (nRet == CE_Failure)
      // return RTN_ERR_GRIDCREATE;
      //
      // GDALClose(pDataSet);
      // delete[] pdRaster;
      // // DEBUG CODE ===========================================================================================================
   }
 
   // // DEBUG CODE ===========================================================================================================
   // string strOutFile = m_strOutPath;
   // strOutFile += "sea_wave_height_before_";
   // strOutFile += to_string(m_ulIter);
   // strOutFile += ".tif";
   //
   // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
   // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
   // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
   // pDataSet->SetGeoTransform(m_dGeoTransform);
   //
   // int nn = 0;
   // double* pdRaster = new double[m_nXGridSize * m_nYGridSize];
   // for (int nY = 0; nY < m_nYGridSize; nY++)
   // {
   // for (int nX = 0; nX < m_nXGridSize; nX++)
   // {
   // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
   // }
   // }
   //
   // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
   // pBand->SetNoDataValue(m_dMissingValue);
   // int nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
   //
   // if (nRet == CE_Failure)
   // return RTN_ERR_GRIDCREATE;
   //
   // GDALClose(pDataSet);
   // delete[] pdRaster;
   // // DEBUG CODE ===========================================================================================================
 
   // // DEBUG CODE ===========================================================================================================
   // strOutFile = m_strOutPath;
   // strOutFile += "sea_wave_angle_before_";
   // strOutFile += to_string(m_ulIter);
   // strOutFile += ".tif";
   //
   // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
   // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
   // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
   // pDataSet->SetGeoTransform(m_dGeoTransform);
   //
   // nn = 0;
   // pdRaster = new double[m_nXGridSize * m_nYGridSize];
   // for (int nY = 0; nY < m_nYGridSize; nY++)
   // {
   // for (int nX = 0; nX < m_nXGridSize; nX++)
   // {
   // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
   // }
   // }
   //
   // pBand = pDataSet->GetRasterBand(1);
   // pBand->SetNoDataValue(m_dMissingValue);
   // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
   //
   // if (nRet == CE_Failure)
   // return RTN_ERR_GRIDCREATE;
   //
   // GDALClose(pDataSet);
   // delete[] pdRaster;
   // // DEBUG CODE ===========================================================================================================
 
   // Now put the X and Y directions together and update the raster cells
   int n = 0;
 
   for (int nY = 0; nY < nYSize; nY++)
   {
      for (int nX = 0; nX < nXSize; nX++)
      {
         int const nActualX = nX + m_nXMinBoundingBox;
         int const nActualY = nY + m_nYMinBoundingBox;
 
         if (m_pRasterGrid->m_Cell[nActualX][nActualY].bIsInContiguousSea())
         {
            // Only update sea cells
            if (m_pRasterGrid->m_Cell[nActualX][nActualY].nGetPolygonID() == INT_NODATA)
            {
               // This is a deep water sea cell (not in a polygon)
               double const dDeepWaterWaveHeight = m_pRasterGrid->m_Cell[nActualX][nActualY].dGetCellDeepWaterWaveHeight();
               m_pRasterGrid->m_Cell[nActualX][nActualY].SetWaveHeight(dDeepWaterWaveHeight);
 
               double const dDeepWaterWaveAngle = m_pRasterGrid->m_Cell[nActualX][nActualY].dGetCellDeepWaterWaveAngle();
               m_pRasterGrid->m_Cell[nActualX][nActualY].SetWaveAngle(dDeepWaterWaveAngle);
            }
 
            else
            {
               // This is in a polygon so is not a deep water sea cell
               double dWaveHeightX;
               double dWaveHeightY;
 
               // Safety checks
               if ((isnan(VdOutX[n])) || (bFPIsEqual(VdOutX[n], m_dMissingValue, TOLERANCE)))
                  dWaveHeightX = dXAvg;
 
               else
                  dWaveHeightX = VdOutX[n];
 
               if ((isnan(VdOutY[n])) || (bFPIsEqual(VdOutY[n], m_dMissingValue, TOLERANCE)))
                  dWaveHeightY = dYAvg;
 
               else
                  dWaveHeightY = VdOutY[n];
 
               // Now calculate wave direction
               double const dWaveHeight = sqrt((dWaveHeightX * dWaveHeightX) + (dWaveHeightY * dWaveHeightY));
               double const dWaveDir = atan2(dWaveHeightX, dWaveHeightY) * (180 / PI);
 
               // assert(isfinite(dWaveHeight));
               // assert(isfinite(dWaveDir));
 
               // Update the cell's wave attributes
               m_pRasterGrid->m_Cell[nActualX][nActualY].SetWaveHeight(dWaveHeight);
               m_pRasterGrid->m_Cell[nActualX][nActualY].SetWaveAngle(dKeepWithin360(dWaveDir));
 
               // Calculate the wave height-to-depth ratio for this cell, then update the cell's active zone status
               double const dSeaDepth = m_pRasterGrid->m_Cell[nActualX][nActualY].dGetSeaDepth();
 
               if ((dWaveHeight / dSeaDepth) >= m_dBreakingWaveHeightDepthRatio)
                  m_pRasterGrid->m_Cell[nActualX][nActualY].SetInActiveZone(true);
 
               // LogStream << " nX = " << nX << " nY = " << nY << " [" << nActualX << "][" << nActualY << "] waveheight = " << dWaveHeight << " dWaveDir = " << dWaveDir << " dKeepWithin360(dWaveDir) = " << dKeepWithin360(dWaveDir) << endl;
            }
         }
 
         // Increment with safety check
         n++;
         n = tMin(n, static_cast<int>(VdOutX.size() - 1));
      }
   }
 
   // // DEBUG CODE ===========================================================================================================
   // strOutFile = m_strOutPath;
   // strOutFile += "sea_wave_height_after_";
   // strOutFile += to_string(m_ulIter);
   // strOutFile += ".tif";
   //
   // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
   // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
   // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
   // pDataSet->SetGeoTransform(m_dGeoTransform);
   //
   // nn = 0;
   // pdRaster = new double[m_nXGridSize * m_nYGridSize];
   // for (int nY = 0; nY < m_nYGridSize; nY++)
   // {
   // for (int nX = 0; nX < m_nXGridSize; nX++)
   // {
   // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveHeight();
   // }
   // }
   //
   // pBand = pDataSet->GetRasterBand(1);
   // pBand->SetNoDataValue(m_dMissingValue);
   // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
   //
   // if (nRet == CE_Failure)
   // return RTN_ERR_GRIDCREATE;
   //
   // GDALClose(pDataSet);
   // delete[] pdRaster;
   // // DEBUG CODE ===========================================================================================================
 
   // // DEBUG CODE ===========================================================================================================
   // strOutFile = m_strOutPath;
   // strOutFile += "sea_wave_angle_after_";
   // strOutFile += to_string(m_ulIter);
   // strOutFile += ".tif";
   //
   // pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
   // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
   // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
   // pDataSet->SetGeoTransform(m_dGeoTransform);
   //
   // nn = 0;
   // pdRaster = new double[m_nXGridSize * m_nYGridSize];
   // for (int nY = 0; nY < m_nYGridSize; nY++)
   // {
   // for (int nX = 0; nX < m_nXGridSize; nX++)
   // {
   // pdRaster[nn++] = m_pRasterGrid->m_Cell[nX][nY].dGetWaveAngle();
   // }
   // }
   //
   // pBand = pDataSet->GetRasterBand(1);
   // pBand->SetNoDataValue(m_dMissingValue);
   // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
   //
   // if (nRet == CE_Failure)
   // return RTN_ERR_GRIDCREATE;
   //
   // GDALClose(pDataSet);
   // delete[] pdRaster;
   // // DEBUG CODE ===========================================================================================================
 
   return RTN_OK;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nInterpolateAllDeepWaterWaveValues(void)
{
   // Interpolate deep water height and orientation from multiple user-supplied values
   unsigned int const nUserPoints = static_cast<unsigned int>(m_VdDeepWaterWaveStationX.size());
 
   // Performance optimization: Enable GDAL threading for interpolation
   CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
 
   // Call GDALGridCreate() with the GGA_InverseDistanceToAPower interpolation algorithm. It has following parameters: radius1 is the first radius (X axis if rotation angle is 0) of the search ellipse, set this to zero (the default) to use the whole point array; radius2 is the second radius (Y axis if rotation angle is 0) of the search ellipse, again set this parameter to zero (the default) to use the whole point array; angle is the angle of the search ellipse rotation in degrees (counter clockwise, default 0.0); nodata is the NODATA marker to fill empty points (default 0.0) TODO 086
   GDALGridInverseDistanceToAPowerOptions *pOptions = new GDALGridInverseDistanceToAPowerOptions();
   pOptions->dfAngle = 0;
   pOptions->dfAnisotropyAngle = 0;
   pOptions->dfAnisotropyRatio = 0;
   pOptions->dfPower = 2;      // Reduced from 3 to 2 for faster computation
   pOptions->dfSmoothing = 50; // Reduced from 100 to 50 for faster computation
   pOptions->dfRadius1 = 0;
   pOptions->dfRadius2 = 0;
   pOptions->nMaxPoints = 12; // Limit points for faster computation (was 0 = unlimited)
   pOptions->nMinPoints = 3;  // Minimum points needed for interpolation
   pOptions->dfNoDataValue = m_nMissingValue;
 
   // CPLSetConfigOption("CPL_DEBUG", "ON");
   // CPLSetConfigOption("GDAL_NUM_THREADS", "1");
 
   // OK, now create a gridded version of wave height: first create the GDAL context TODO 086
   // GDALGridContext* pContext = GDALGridContextCreate(GGA_InverseDistanceToAPower, pOptions, nUserPoints, &m_VdDeepWaterWaveStationX[0], &m_VdDeepWaterWaveStationY[0], &m_VdThisIterDeepWaterWaveStationHeight[0], true);
   GDALGridContext *pContext = GDALGridContextCreate(GGA_InverseDistanceToAPower, pOptions, nUserPoints, m_VdDeepWaterWaveStationX.data(), m_VdDeepWaterWaveStationY.data(), m_VdThisIterDeepWaterWaveStationHeight.data(), true);
 
   if (pContext == NULL)
   {
      delete pOptions;
      return RTN_ERR_GRIDCREATE;
   }
 
   // Now process the context
   double *dHeightOut = new double[m_ulNumCells];
   int nRet = GDALGridContextProcess(pContext, 0, m_nXGridSize - 1, 0, m_nYGridSize - 1, m_nXGridSize, m_nYGridSize, GDT_Float64, dHeightOut, NULL, NULL);
 
   if (nRet == CE_Failure)
   {
      delete[] dHeightOut;
      delete pOptions;
      return RTN_ERR_GRIDCREATE;
   }
 
   // Get rid of the context
   GDALGridContextFree(pContext);
 
   // Next create a gridded version of wave orientation: first create the GDAL context
   // pContext = GDALGridContextCreate(GGA_InverseDistanceToAPower, pOptions, nUserPoints,  &(m_VdDeepWaterWaveStationX[0]), &(m_VdDeepWaterWaveStationY[0]), (&m_VdThisIterDeepWaterWaveStationAngle[0]), true);
   pContext = GDALGridContextCreate(GGA_InverseDistanceToAPower, pOptions, nUserPoints, m_VdDeepWaterWaveStationX.data(), m_VdDeepWaterWaveStationY.data(), m_VdThisIterDeepWaterWaveStationAngle.data(), true);
 
   if (pContext == NULL)
   {
      delete[] dHeightOut;
      delete pOptions;
      return RTN_ERR_GRIDCREATE;
   }
 
   // Now process the context TODO 086
   double *dAngleOut = new double[m_ulNumCells];
   nRet = GDALGridContextProcess(pContext, 0, m_nXGridSize - 1, 0, m_nYGridSize - 1, m_nXGridSize, m_nYGridSize, GDT_Float64, dAngleOut, NULL, NULL);
 
   if (nRet == CE_Failure)
   {
      delete[] dHeightOut;
      delete[] dAngleOut;
      delete pOptions;
      return RTN_ERR_GRIDCREATE;
   }
 
   // Get rid of the context
   GDALGridContextFree(pContext);
 
   // OK, now create a gridded version of wave period: first create the GDAL context
   // pContext = GDALGridContextCreate(GGA_InverseDistanceToAPower, pOptions, nUserPoints, &m_VdDeepWaterWaveStationX[0], &m_VdDeepWaterWaveStationY[0], &m_VdThisIterDeepWaterWaveStationPeriod[0], true);
   pContext = GDALGridContextCreate(GGA_InverseDistanceToAPower, pOptions, nUserPoints, m_VdDeepWaterWaveStationX.data(), m_VdDeepWaterWaveStationY.data(), m_VdThisIterDeepWaterWaveStationPeriod.data(), true);
 
   if (pContext == NULL)
   {
      delete pOptions;
      return RTN_ERR_GRIDCREATE;
   }
 
   // Now process the context TODO 086
   double *dPeriopdOut = new double[m_ulNumCells];
   nRet = GDALGridContextProcess(pContext, 0, m_nXGridSize - 1, 0, m_nYGridSize - 1, m_nXGridSize, m_nYGridSize, GDT_Float64, dPeriopdOut, NULL, NULL);
 
   if (nRet == CE_Failure)
   {
      delete[] dPeriopdOut;
      delete pOptions;
      return RTN_ERR_GRIDCREATE;
   }
 
   // Get rid of the context
   GDALGridContextFree(pContext);
 
   // The output from GDALGridCreate() is in dHeightOut, dAngleOut and dPeriopdOut but must be reversed
   vector<double> VdHeight;
   vector<double> VdAngle;
   vector<double> VdPeriod;
 
   int n = 0;
   int nValidHeight = 0;
   int nValidAngle = 0;
   int nValidPeriod = 0;
 
   double dAvgHeight = 0;
   double dAvgAngle = 0;
   double dAvgPeriod = 0;
 
   for (int nY = m_nYGridSize - 1; nY >= 0; nY--)
   {
      for (int nX = 0; nX < m_nXGridSize; nX++)
      {
         if (isfinite(dHeightOut[n]))
         {
            VdHeight.push_back(dHeightOut[n]);
 
            dAvgHeight += dHeightOut[n];
            nValidHeight++;
         }
 
         else
         {
            VdHeight.push_back(m_dMissingValue);
         }
 
         if (isfinite(dAngleOut[n]))
         {
            VdAngle.push_back(dAngleOut[n]);
 
            dAvgAngle += dAngleOut[n];
            nValidAngle++;
         }
 
         else
         {
            VdAngle.push_back(m_dMissingValue);
         }
 
         if (isfinite(dPeriopdOut[n]))
         {
            VdPeriod.push_back(dPeriopdOut[n]);
 
            dAvgPeriod += dPeriopdOut[n];
            nValidPeriod++;
         }
 
         else
         {
            VdPeriod.push_back(m_dMissingValue);
         }
 
         // LogStream << " nX = " << nX << " nY = " << nY << " n = " << n << " dHeightOut[n] = " << dHeightOut[n] << " dAngleOut[n] = " << dAngleOut[n] << endl;
         n++;
      }
   }
 
   // Calculate averages
   dAvgHeight /= nValidHeight;
   dAvgAngle /= nValidAngle;
   dAvgPeriod /= nValidPeriod;
 
   // Tidy
   delete pOptions;
   delete[] dHeightOut;
   delete[] dAngleOut;
   delete[] dPeriopdOut;
 
   // Now update all raster cells
   n = 0;
 
   for (int nY = 0; nY < m_nYGridSize; nY++)
   {
      for (int nX = 0; nX < m_nXGridSize; nX++)
      {
         if (bFPIsEqual(VdHeight[n], m_dMissingValue, TOLERANCE))
            m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWaveHeight(dAvgHeight);
 
         else
            m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWaveHeight(VdHeight[n]);
 
         if (bFPIsEqual(VdAngle[n], m_dMissingValue, TOLERANCE))
            m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWaveAngle(dAvgAngle);
 
         else
            m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWaveAngle(VdAngle[n]);
 
         if (bFPIsEqual(VdPeriod[n], m_dMissingValue, TOLERANCE))
            m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWavePeriod(dAvgPeriod);
 
         else
            m_pRasterGrid->m_Cell[nX][nY].SetCellDeepWaterWavePeriod(VdPeriod[n]);
 
         // LogStream << " [" << nX << "][" << nY << "] deep water wave height = " << m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveHeight() << " deep water wave angle = " << m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveAngle() << endl;
         n++;
      }
   }
 
   // // DEBUG CODE ===========================================================================================================
   // string strOutFile = m_strOutPath;
   // strOutFile += "init_deep_water_wave_height_";
   // strOutFile += to_string(m_ulIter);
   // strOutFile += ".tif";
   // GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("gtiff");
   // GDALDataset* pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
   // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
   // pDataSet->SetGeoTransform(m_dGeoTransform);
   // double* pdRaster = new double[m_ulNumCells];
   // int nn = 0;
   // for (int nY = 0; nY < m_nYGridSize; nY++)
   // {
   // for (int nX = 0; nX < m_nXGridSize; nX++)
   // {
   //          // Write this value to the array
   // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveHeight();
   // nn++;
   // }
   // }
   //
   // GDALRasterBand* pBand = pDataSet->GetRasterBand(1);
   // pBand->SetNoDataValue(m_nMissingValue);
   // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
   //
   // if (nRet == CE_Failure)
   // return RTN_ERR_GRIDCREATE;
   //
   // GDALClose(pDataSet);
   // // DEBUG CODE ===========================================================================================================
 
   // // DEBUG CODE ===========================================================================================================
   // strOutFile = m_strOutPath;
   // strOutFile += "init_deep_water_wave_angle_";
   // strOutFile += to_string(m_ulIter);
   // strOutFile += ".tif";
   // pDataSet = pDriver->Create(strOutFile.c_str(), m_nXGridSize, m_nYGridSize, 1, GDT_Float64, m_papszGDALRasterOptions);
   // pDataSet->SetProjection(m_strGDALBasementDEMProjection.c_str());
   // pDataSet->SetGeoTransform(m_dGeoTransform);
   // nn = 0;
   // for (int nY = 0; nY < m_nYGridSize; nY++)
   // {
   // for (int nX = 0; nX < m_nXGridSize; nX++)
   // {
   //          // Write this value to the array
   // pdRaster[nn] = m_pRasterGrid->m_Cell[nX][nY].dGetCellDeepWaterWaveAngle();
   // nn++;
   // }
   // }
   //
   // pBand = pDataSet->GetRasterBand(1);
   // pBand->SetNoDataValue(m_nMissingValue);
   // nRet = pBand->RasterIO(GF_Write, 0, 0, m_nXGridSize, m_nYGridSize, pdRaster, m_nXGridSize, m_nYGridSize, GDT_Float64, 0, 0, NULL);
   //
   // if (nRet == CE_Failure)
   // return RTN_ERR_GRIDCREATE;
   //
   // GDALClose(pDataSet);
   // delete[] pdRaster;
   // // DEBUG CODE ===========================================================================================================
 
   return RTN_OK;
}