1 /*!
2 *
3 * \file gis_utils.cpp
4 * \brief Various GIS-related functions. This version will build with GDAL version 2
5 * \details TODO A more detailed description of these routines.
6 * \author Andres Payo
7 * \author David Favis-Mortlock
8 * \author Martin Husrt
9 * \author Monica Palaseanu-Lovejoy
10 * \date 2017
11 * \copyright GNU General Public License
12 *
13 */
14
15 /*===============================================================================================================================
16
17 This file is part of CliffMetrics, the Coastal Modelling Environment.
18
19 CliffMetrics 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.
20
21 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.
22
23 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.
24
25 ===============================================================================================================================*/
26 #include <iostream>
27 using std::cout;
28 using std::cerr;
29 using std::endl;
30 using std::ios;
31
32 #include <cmath>
33 #include <cfloat>
34
35 #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
36 #include <gdal_priv.h>
37 #include <ogrsf_frmts.h>
38 #endif // #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
39
40 #include "cliffmetrics.h"
41 #include "delineation.h"
42 #include "raster_grid.h"
43
44
45 /*
46 Notes re. co-ordinate systems used
47
48 1. In the raster CRS, cell[0][0] is at the top left (NW) corner of the grid. Raster grid co-oordinate [0][0] is actually the top left (NW) corner of this cell.
49
50 2. We assume that the grid CRS and external CRS have parallel axes. If they have not, see http://www.gdal.org/classGDALDataset.html which says that:
51
52 To convert between pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space
53 Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
54 Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];
55
56 In a north-up image, padfTransform[1] is the pixel width, and padfTransform[5] is the pixel height. The upper left corner of the upper left pixel is at position
57 (padfTransform[0], padfTransform[3]).
58
59 3. Usually, raster grid CRS values are integer, i.e. they refer to a point which is at the centroid of a cell. They may also be -ve or greater than m_nXGridMax-1 i.e. may refer to a point which lies outside any cell of the raster grid.
60
61 */
62
63 /*==============================================================================================================================
64
65 Given the integer X-axis ordinate of a cell in the raster-grid CRS, returns the external CRS X-axis ordinate of the cell's centroid
66
67 ===============================================================================================================================*/
dGridCentroidXToExtCRSX(int const nGridX) const68 double CDelineation::dGridCentroidXToExtCRSX(int const nGridX) const
69 {
70 return (m_dGeoTransform[0] + (nGridX * m_dGeoTransform[1]) + (m_dGeoTransform[1] / 2));
71 }
72
73
74 /*==============================================================================================================================
75
76 Given the integer Y-axis ordinate of a cell in the raster-grid CRS, returns the external CRS Y-axis ordinate of the cell's centroid
77
78 ===============================================================================================================================*/
dGridCentroidYToExtCRSY(int const nGridY) const79 double CDelineation::dGridCentroidYToExtCRSY(int const nGridY) const
80 {
81 return (m_dGeoTransform[3] + (nGridY * m_dGeoTransform[5]) + (m_dGeoTransform[5] / 2));
82 }
83
84 /*==============================================================================================================================
85
86 Transforms a pointer to a C2DIPoint in the raster-grid CRS (assumed to be the centroid of a cell) to the equivalent C2DPoint in the external CRS
87
88 ===============================================================================================================================*/
PtGridCentroidToExt(C2DIPoint * const pPtiIn) const89 C2DPoint CDelineation::PtGridCentroidToExt(C2DIPoint* const pPtiIn) const
90 {
91 int
92 nGridX = pPtiIn->nGetX(),
93 nGridY = pPtiIn->nGetY();
94
95 double
96 dX = m_dGeoTransform[0] + (nGridX * m_dGeoTransform[1]) + (m_dGeoTransform[1] / 2),
97 dY = m_dGeoTransform[3] + (nGridY * m_dGeoTransform[5]) + (m_dGeoTransform[5] / 2);
98
99 return C2DPoint(dX, dY);
100 }
101
102
103 /*==============================================================================================================================
104
105 Given a real-valued X-axis ordinate in the raster-grid CRS (i.e. not the centroid of a cell), returns the external CRS X-axis ordinate
106
107 ===============================================================================================================================*/
dGridXToExtCRSX(double const dGridX) const108 double CDelineation::dGridXToExtCRSX(double const dGridX) const
109 {
110 return (m_dGeoTransform[0] + (dGridX * m_dGeoTransform[1]));
111 }
112
113
114 /*==============================================================================================================================
115
116 Given a real-valued Y-axis ordinate in the raster-grid CRS (i.e. not the centroid of a cell), returns the external CRS Y-axis ordinate
117
118 ===============================================================================================================================*/
dGridYToExtCRSY(double const dGridY) const119 double CDelineation::dGridYToExtCRSY(double const dGridY) const
120 {
121 return (m_dGeoTransform[3] + (dGridY * m_dGeoTransform[5]));
122 }
123
124
125 /*==============================================================================================================================
126
127 Transforms an X-axis ordinate in the external CRS to the equivalent X-axis ordinate in the raster-grid CRS (the result may not be integer, and may be outside the grid)
128
129 ===============================================================================================================================*/
dExtCRSXToGridX(double const dExtCRSX) const130 double CDelineation::dExtCRSXToGridX(double const dExtCRSX) const
131 {
132 return ((dExtCRSX - m_dGeoTransform[0]) / m_dGeoTransform[1]);
133 }
134
135
136 /*==============================================================================================================================
137
138 Transforms a Y-axis ordinate in the external CRS to the equivalent Y-axis ordinate in the raster-grid CRS (the result may not be integer, and may be outside the grid)
139
140 ===============================================================================================================================*/
dExtCRSYToGridY(double const dExtCRSY) const141 double CDelineation::dExtCRSYToGridY(double const dExtCRSY) const
142 {
143 return ((dExtCRSY - m_dGeoTransform[3]) / m_dGeoTransform[5]);
144 }
145
146
147 /*==============================================================================================================================
148
149 Transforms a pointer to a C2DPoint in the external CRS to the equivalent C2DIPoint in the raster-grid CRS (both values rounded)
150
151 ===============================================================================================================================*/
PtiExtCRSToGrid(C2DPoint * const pPtIn) const152 C2DIPoint CDelineation::PtiExtCRSToGrid(C2DPoint* const pPtIn) const
153 {
154 double
155 dX = pPtIn->dGetX(),
156 dY = pPtIn->dGetY();
157
158 int
159 nX = dRound((dX - m_dGeoTransform[0]) / m_dGeoTransform[1]),
160 nY = dRound((dY - m_dGeoTransform[3]) / m_dGeoTransform[5]);
161
162 return C2DIPoint(nX, nY);
163 }
164
165
166 /*==============================================================================================================================
167
168 Returns the distance (in external CRS) between two points
169
170 ===============================================================================================================================*/
dGetDistanceBetween(C2DPoint * const Pt1,C2DPoint * const Pt2)171 double CDelineation::dGetDistanceBetween(C2DPoint* const Pt1, C2DPoint* const Pt2)
172 {
173 double
174 dXDist = Pt1->dGetX() - Pt2->dGetX(),
175 dYDist = Pt1->dGetY() - Pt2->dGetY();
176
177 return hypot(dXDist, dYDist);
178 }
179
180
181 /*==============================================================================================================================
182
183 Returns the distance (in grid units) between two grid cell points
184
185 ===============================================================================================================================*/
dGetDistanceBetween(C2DIPoint * const Pti1,C2DIPoint * const Pti2)186 double CDelineation::dGetDistanceBetween(C2DIPoint* const Pti1, C2DIPoint* const Pti2)
187 {
188 double
189 dXDist = Pti1->nGetX() - Pti2->nGetX(),
190 dYDist = Pti1->nGetY() - Pti2->nGetY();
191
192 return hypot(dXDist, dYDist);
193 }
194
195
196 /*==============================================================================================================================
197
198 Checks whether the supplied point (an x-y pair, in the grid CRS) is within the raster grid
199
200 ===============================================================================================================================*/
bIsWithinGrid(int const nX,int const nY) const201 bool CDelineation::bIsWithinGrid(int const nX, int const nY) const
202 {
203 if ((nX < 0) || (nX >= m_nXGridMax) || (nY < 0) || (nY >= m_nYGridMax))
204 return false;
205
206 return true;
207 }
208
209
210 /*==============================================================================================================================
211
212 Checks whether the supplied point (a reference to a C2DIPoint, in the grid CRS) is within the raster grid
213
214 ===============================================================================================================================*/
bIsWithinGrid(C2DIPoint * const Pti) const215 bool CDelineation::bIsWithinGrid(C2DIPoint* const Pti) const
216 {
217 int nX = Pti->nGetX();
218
219 if ((nX < 0) || (nX >= m_nXGridMax))
220 return false;
221
222 int nY = Pti->nGetY();
223
224 if ((nY < 0) || (nY >= m_nYGridMax))
225 return false;
226
227 return true;
228 }
229
230
231 /*==============================================================================================================================
232
233 Constrains the supplied point (a reference to a C2DIPoint, in the grid CRS) to be within the raster grid
234
235 ===============================================================================================================================*/
KeepWithinGrid(C2DIPoint * Pti)236 void CDelineation::KeepWithinGrid(C2DIPoint* Pti)
237 {
238 int nX = Pti->nGetX();
239 nX = tMax(nX, 0);
240 nX = tMin(nX, m_nXGridMax-1);
241 Pti->SetX(nX);
242
243 int nY = Pti->nGetY();
244 nY = tMax(nY, 0);
245 nY = tMin(nY, m_nYGridMax-1);
246 Pti->SetY(nY);
247 }
248
249
250 /*==============================================================================================================================
251
252 Constrains the supplied point (an x-y pair, in the grid CRS) to be within the raster grid
253
254 ===============================================================================================================================*/
KeepWithinGrid(int & nX,int & nY)255 void CDelineation::KeepWithinGrid(int& nX, int& nY)
256 {
257 nX = tMax(nX, 0);
258 nX = tMin(nX, m_nXGridMax-1);
259
260 nY = tMax(nY, 0);
261 nY = tMin(nY, m_nYGridMax-1);
262 }
263
264
265 /*==============================================================================================================================
266
267 Constrains the supplied point (a reference to a C2DPoint, in the external CRS) to be within the raster grid
268
269 ===============================================================================================================================*/
270 // C2DPoint* CDelineation::pPtExtCRSKeepWithinGrid(C2DPoint* pPt)
271 // {
272 // double dGridX = dExtCRSXToGridX(pPt->dGetX());
273 // dGridX = tMax(dGridX, 0.0);
274 // dGridX = tMin(dGridX, static_cast<double>(m_nXGridMax-1));
275 // pPt->SetX(dGridXToExtCRSX(dGridX));
276 //
277 // double dGridY = dExtCRSYToGridY(pPt->dGetY());
278 // dGridY = tMax(dGridY, 0.0);
279 // dGridY = tMin(dGridY, static_cast<double>(m_nYGridMax-1));
280 // pPt->SetY(dGridYToExtCRSY(dGridY));
281 //
282 // return pPt;
283 // }
284
285
286 /*==============================================================================================================================
287
288 Constrains the supplied angle to be within 0 and 360 degrees
289
290 ===============================================================================================================================*/
dKeepWithin360(double const dAngle)291 double CDelineation::dKeepWithin360(double const dAngle)
292 {
293 double dNewAngle = fmod(dAngle, 360);
294
295 // Sort out -ve angles
296 while (dNewAngle < 0)
297 dNewAngle += 360;
298
299 return dNewAngle;
300 }
301
302
303 /*==============================================================================================================================
304
305 Returns a point (external CRS) which is the average of (i.e. is midway between) two other external CRS points
306
307 ===============================================================================================================================*/
PtAverage(C2DPoint * const pPt1,C2DPoint * const pPt2)308 C2DPoint CDelineation::PtAverage(C2DPoint* const pPt1, C2DPoint* const pPt2)
309 {
310 double
311 dPt1X = pPt1->dGetX(),
312 dPt1Y = pPt1->dGetY(),
313 dPt2X = pPt2->dGetX(),
314 dPt2Y = pPt2->dGetY(),
315 dPtAvgX = (dPt1X + dPt2X) / 2,
316 dPtAvgY = (dPt1Y + dPt2Y) / 2;
317
318 return C2DPoint(dPtAvgX, dPtAvgY);
319 }
320
321
322 /*==============================================================================================================================
323
324 Returns a point (external CRS) which is the average of a vector of external CRS points
325
326 ===============================================================================================================================*/
PtAverage(vector<C2DPoint> * pVIn)327 C2DPoint CDelineation::PtAverage(vector<C2DPoint>* pVIn)
328 {
329 int nSize = pVIn->size();
330 if (nSize == 0)
331 return C2DPoint(DBL_NODATA, DBL_NODATA);
332
333 double
334 dAvgX = 0,
335 dAvgY = 0;
336
337 for (int n = 0; n < nSize; n++)
338 {
339 dAvgX += pVIn->at(n).dGetX();
340 dAvgY += pVIn->at(n).dGetY();
341 }
342
343 dAvgX /= nSize;
344 dAvgY /= nSize;
345
346 return C2DPoint(dAvgX, dAvgY);
347 }
348
349
350 /*==============================================================================================================================
351
352 Returns a point (grid CRS) which is the rounded average of (i.e. is midway between) two other grid CRS points
353
354 ===============================================================================================================================*/
PtiAverage(C2DIPoint * const pPti1,C2DIPoint * const pPti2)355 C2DIPoint CDelineation::PtiAverage(C2DIPoint* const pPti1, C2DIPoint* const pPti2)
356 {
357 double
358 dPt1X = pPti1->nGetX(),
359 dPt1Y = pPti1->nGetY(),
360 dPt2X = pPti2->nGetX(),
361 dPt2Y = pPti2->nGetY(),
362 dPtAvgX = (dPt1X + dPt2X) / 2,
363 dPtAvgY = (dPt1Y + dPt2Y) / 2;
364
365 return C2DIPoint(static_cast<int>(dRound(dPtAvgX)), static_cast<int>(dRound(dPtAvgY)));
366 }
367
368
369 /*==============================================================================================================================
370
371 Returns a vector which is perpendicular to an existing vector
372
373 ===============================================================================================================================*/
374 // vector<C2DPoint> CDelineation::VGetPerpendicular(C2DPoint* const PtStart, C2DPoint* const PtNext, double const dDesiredLength, int const nHandedness)
375 // {
376 // // Returns a two-point vector which passes through PtStart with a scaled length
377 // double dXLen = PtNext->dGetX() - PtStart->dGetX();
378 // double dYLen = PtNext->dGetY() - PtStart->dGetY();
379 //
380 // double dLength = hypot(dXLen, dYLen);
381 // double dScaleFactor = dDesiredLength / dLength;
382 //
383 // // The difference vector is (dXLen, dYLen), so the perpendicular difference vector is (-dYLen, dXLen) or (dYLen, -dXLen)
384 // C2DPoint EndPt;
385 // if (nHandedness == RIGHT_HANDED)
386 // {
387 // EndPt.SetX(PtStart->dGetX() + (dScaleFactor * dYLen));
388 // EndPt.SetY(PtStart->dGetY() - (dScaleFactor * dXLen));
389 // }
390 // else
391 // {
392 // EndPt.SetX(PtStart->dGetX() - (dScaleFactor * dYLen));
393 // EndPt.SetY(PtStart->dGetY() + (dScaleFactor * dXLen));
394 // }
395 //
396 // vector<C2DPoint> VNew;
397 // VNew.push_back(*PtStart);
398 // VNew.push_back(EndPt);
399 // return VNew;
400 // }
401
402
403
404 /*==============================================================================================================================
405
406 Returns a C2DPoint which is the 'other' point of a two-point vector passing through PtStart, and which is perpendicular to the two-point vector from PtStart to PtNext
407
408 ===============================================================================================================================*/
PtGetPerpendicular(C2DPoint * const PtStart,C2DPoint * const PtNext,double const dDesiredLength,int const nHandedness)409 C2DPoint CDelineation::PtGetPerpendicular(C2DPoint* const PtStart, C2DPoint* const PtNext, double const dDesiredLength, int const nHandedness)
410 {
411 double dXLen = PtNext->dGetX() - PtStart->dGetX();
412 double dYLen = PtNext->dGetY() - PtStart->dGetY();
413
414 double dLength = hypot(dXLen, dYLen);
415 double dScaleFactor = dDesiredLength / dLength;
416
417 // The difference vector is (dXLen, dYLen), so the perpendicular difference vector is (-dYLen, dXLen) or (dYLen, -dXLen)
418 C2DPoint EndPt;
419 if (nHandedness == RIGHT_HANDED)
420 {
421 EndPt.SetX(PtStart->dGetX() + (dScaleFactor * dYLen));
422 EndPt.SetY(PtStart->dGetY() - (dScaleFactor * dXLen));
423 }
424 else
425 {
426 EndPt.SetX(PtStart->dGetX() - (dScaleFactor * dYLen));
427 EndPt.SetY(PtStart->dGetY() + (dScaleFactor * dXLen));
428 }
429
430 return EndPt;
431 }
432
433
434 /*==============================================================================================================================
435
436 Checks whether the selected raster GDAL driver supports file creation, 32-bit doubles, etc.
437
438 ===============================================================================================================================*/
bCheckRasterGISOutputFormat(void)439 bool CDelineation::bCheckRasterGISOutputFormat(void)
440 {
441 #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
442 // Register all available GDAL raster and vector drivers (GDAL 2)
443 GDALAllRegister();
444
445 // If the user hasn't specified a GIS output format, assume that we will use the same GIS format as the input basement DEM
446 if (m_strRasterGISOutFormat.empty())
447 m_strRasterGISOutFormat = m_strGDALBasementDEMDriverCode;
448
449 // Load the raster GDAL driver
450 GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName(m_strRasterGISOutFormat.c_str());
451 if (NULL == pDriver)
452 {
453 // Can't load raster GDAL driver. Incorrectly specified?
454 cerr << ERR << "Unknown raster GIS output format '" << m_strRasterGISOutFormat << "'." << endl;
455 return false;
456 }
457
458 // Get the metadata for this raster driver
459 char** papszMetadata = pDriver->GetMetadata();
460
461 // for (int i = 0; papszMetadata[i] != NULL; i++)
462 // cout << papszMetadata[i] << endl;
463 // cout << endl;
464
465 // For GDAL2, need to test if this is a raster driver
466 if (! CSLFetchBoolean(papszMetadata, GDAL_DCAP_RASTER, FALSE))
467 {
468 // This is not a raster driver
469 cerr << ERR << "GDAL driver '" << m_strRasterGISOutFormat << "' is not a raster driver. Choose another format." << endl;
470 return false;
471 }
472
473 // This driver is OK, so store its longname and the default file extension
474 m_strGDALRasterOutputDriverLongname = CSLFetchNameValue(papszMetadata, "DMD_LONGNAME");
475 m_strGDALRasterOutputDriverExtension = CSLFetchNameValue(papszMetadata, "DMD_EXTENSION");
476
477 // Set up any defaults for raster files that are created using this driver
478 SetRasterFileCreationDefaults();
479
480 // Now do various tests of the driver's capabilities
481 if (! CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
482 {
483 // This raster driver does not support the Create() method, does it support CreateCopy()?
484 if (! CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATECOPY, FALSE))
485 {
486 cerr << ERR << "Cannot write using raster GDAL driver '" << m_strRasterGISOutFormat << " since neither Create() or CreateCopy() are supported'. Choose another GDAL raster format." << endl;
487 return false;
488 }
489
490 // Can't use Create() but can use CreateCopy()
491 m_bGDALCanCreate = false;
492 }
493
494 // Next, test to see what data types the driver can write and from this, work out the largest int and float we can write
495 if (strstr(CSLFetchNameValue(papszMetadata, "DMD_CREATIONDATATYPES"), "Float"))
496 {
497 m_bGDALCanWriteFloat = true;
498 m_GDALWriteFloatDataType = GDT_Float32;
499 }
500
501 if (strstr(CSLFetchNameValue(papszMetadata, "DMD_CREATIONDATATYPES"), "UInt32"))
502 {
503 m_bGDALCanWriteInt32 = true;
504
505 m_GDALWriteIntDataType = GDT_UInt32;
506 m_lGDALMaxCanWrite = UINT32_MAX;
507 m_lGDALMinCanWrite = 0;
508
509 if (! m_bGDALCanWriteFloat)
510 m_GDALWriteFloatDataType = GDT_UInt32;
511
512 return true;
513 }
514
515 if (strstr(CSLFetchNameValue(papszMetadata, "DMD_CREATIONDATATYPES"), "Int32"))
516 {
517 m_bGDALCanWriteInt32 = true;
518
519 m_GDALWriteIntDataType = GDT_Int32;
520 m_lGDALMaxCanWrite = INT32_MAX;
521 m_lGDALMinCanWrite = INT32_MIN;
522
523 if (! m_bGDALCanWriteFloat)
524 m_GDALWriteFloatDataType = GDT_Int32;
525
526 return true;
527 }
528
529 if (strstr(CSLFetchNameValue(papszMetadata, "DMD_CREATIONDATATYPES"), "UInt16"))
530 {
531 m_bGDALCanWriteInt32 = false;
532
533 m_GDALWriteIntDataType = GDT_UInt16;
534 m_lGDALMaxCanWrite = UINT16_MAX;
535 m_lGDALMinCanWrite = 0;
536
537 if (! m_bGDALCanWriteFloat)
538 m_GDALWriteFloatDataType = GDT_UInt16;
539
540 return true;
541 }
542
543 if (strstr(CSLFetchNameValue(papszMetadata, "DMD_CREATIONDATATYPES"), "Int16"))
544 {
545 m_bGDALCanWriteInt32 = false;
546
547 m_GDALWriteIntDataType = GDT_Int16;
548 m_lGDALMaxCanWrite = INT16_MAX;
549 m_lGDALMinCanWrite = INT16_MIN;
550
551 if (! m_bGDALCanWriteFloat)
552 m_GDALWriteFloatDataType = GDT_Int16;
553
554 return true;
555 }
556
557 if (strstr(CSLFetchNameValue(papszMetadata, "DMD_CREATIONDATATYPES"), "Byte"))
558 {
559 m_bGDALCanWriteInt32 = false;
560
561 m_GDALWriteIntDataType = GDT_Byte;
562 m_lGDALMaxCanWrite = UINT8_MAX;
563 m_lGDALMinCanWrite = 0;
564
565 if (! m_bGDALCanWriteFloat)
566 m_GDALWriteFloatDataType = GDT_Byte;
567
568 return true;
569 }
570
571 // This driver does not even support byte output
572 cerr << ERR << "Cannot write using raster GDAL driver '" << m_strRasterGISOutFormat << ", not even byte output is supported'. Choose another GIS raster format." << endl;
573 return false;
574 #else // #if defined(_SAGA_MSW) || defined(_SAGA_LINUX)
575 return true;
576 #endif // #if defined(_SAGA_MSW) || defined(_SAGA_LINUX)
577 }
578
579
580 /*==============================================================================================================================
581
582 Checks whether the selected vector OGR driver supports file creation etc.
583
584 ===============================================================================================================================*/
bCheckVectorGISOutputFormat(void)585 bool CDelineation::bCheckVectorGISOutputFormat(void)
586 {
587 #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
588 // Load the vector GDAL driver (NOTE this assumes that GDALAllRegister() has already been called)
589 GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName(m_strVectorGISOutFormat.c_str());
590 if (NULL == pDriver)
591 {
592 // Can't load vector GDAL driver. Incorrectly specified?
593 cerr << ERR << "Unknown vector GIS output format '" << m_strVectorGISOutFormat << "'." << endl;
594 return false;
595 }
596
597 // Get the metadata for this vector driver
598 char** papszMetadata = pDriver->GetMetadata();
599
600 // For GDAL2, need to test if this is a vector driver
601 if (! CSLFetchBoolean(papszMetadata, GDAL_DCAP_VECTOR, FALSE))
602 {
603 // This is not a vector driver
604 cerr << ERR << "GDAL driver '" << m_strVectorGISOutFormat << "' is not a vector driver. Choose another format." << endl;
605 return false;
606 }
607
608 if (! CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
609 {
610 // Driver does not support create() method
611 cerr << ERR << "Cannot write vector GIS files using GDAL driver '" << m_strRasterGISOutFormat << "'. Choose another format." << endl;
612 return false;
613 }
614
615 // Driver is OK, now set some options for individual drivers
616 if (m_strVectorGISOutFormat == "ESRI Shapefile")
617 {
618 // Set this, so that just a single dataset-with-one-layer shapefile is created, rather than a directory
619 // (see http://www.gdal.org/ogr/drv_shapefile.html)
620 m_strOGRVectorOutputExtension = ".shp";
621
622 }
623 // TODO Others
624 #endif // #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
625
626 return true;
627 }
628
629
630 /*==============================================================================================================================
631
632 The bSaveAllRasterGISFiles member function saves the raster GIS files using values from the RasterGrid array
633
634 ==============================================================================================================================*/
bSaveAllRasterGISFiles(void)635 bool CDelineation::bSaveAllRasterGISFiles(void)
636 {
637 #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
638 // Increment file number
639 // m_nGISSave++;
640
641 // Set for next save
642 // m_nThisSave = tMin(++m_nThisSave, m_nUSave);
643
644 // These are always written
645 if (! bWriteRasterGISFloat(PLOT_SEDIMENT_TOP_ELEV, &PLOT_SEDIMENT_TOP_ELEV_TITLE))
646 return false;
647
648 if (m_bRasterCoastlineSave)
649 {
650 if (! bWriteRasterGISInt(PLOT_RASTER_COAST, &PLOT_RASTER_COAST_TITLE))
651 return false;
652 }
653
654 if (m_bRasterNormalSave)
655 {
656 if (! bWriteRasterGISInt(PLOT_RASTER_NORMAL, &PLOT_RASTER_NORMAL_TITLE))
657 return false;
658 }
659
660 // These are optional
661 //if (! bWriteRasterGISInt(PLOT_LANDFORM, &PLOT_LANDFORM_TITLE))
662 // return false;*/
663
664 return true;
665 #else // #if defined(_SAGA_MSW) || defined(_SAGA_LINUX)
666 CSG_Grid *pGrid;
667
668 if( (pGrid = (*m_pParameters)("SEDIMENT_TOP" )->asGrid()) )
669 {
670 if( !bWriteRasterGISFloat(PLOT_SEDIMENT_TOP_ELEV, pGrid) )
671 {
672 return( false );
673 }
674 }
675
676 if( (pGrid = (*m_pParameters)("RASTER_COAST" )->asGrid()) )
677 {
678 if( !bWriteRasterGISInt (PLOT_RASTER_COAST , pGrid) )
679 {
680 return( false );
681 }
682 }
683
684 if( (pGrid = (*m_pParameters)("RASTER_NORMAL")->asGrid()) )
685 {
686 if( !bWriteRasterGISInt (PLOT_RASTER_NORMAL , pGrid) )
687 {
688 return( false );
689 }
690 }
691
692 //if( (pGrid = (*m_pParameters)("PLOT_LANDFORM")->asGrid()) )
693 //{
694 // if( !bWriteRasterGISFloat(PLOT_LANDFORM , pGrid) )
695 // {
696 // return( false );
697 // }
698 //}
699
700 return( true );
701 #endif // #if defined(_SAGA_MSW) || defined(_SAGA_LINUX)
702 }
703
704
705 /*==============================================================================================================================
706
707 The bSaveAllvectorGISFiles member function saves the vector GIS files
708
709 ==============================================================================================================================*/
bSaveAllVectorGISFiles(void)710 bool CDelineation::bSaveAllVectorGISFiles(void)
711 {
712 #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
713 // Always written
714 if (! bWriteVectorGIS(PLOT_COAST, &PLOT_COAST_TITLE))
715 return false;
716
717 if (! bWriteVectorGIS(PLOT_NORMALS, &PLOT_NORMALS_TITLE))
718 return false;
719
720 if (! bWriteVectorGIS(PLOT_CLIFF_TOP, &PLOT_CLIFF_TOP_TITLE))
721 return false;
722
723 if (! bWriteVectorGIS(PLOT_CLIFF_TOE, &PLOT_CLIFF_TOE_TITLE))
724 return false;
725
726 if (! bWriteVectorGIS(PLOT_COAST_POINT, &PLOT_COAST_POINT_TITLE))
727 return false;
728
729 // These are optional (see delineation constructor)
730 if (m_bInvalidNormalsSave)
731 {
732 if (! bWriteVectorGIS(PLOT_INVALID_NORMALS, &PLOT_INVALID_NORMALS_TITLE))
733 return false;
734 }
735
736 if (m_bCoastCurvatureSave)
737 {
738 if (! bWriteVectorGIS(PLOT_COAST_CURVATURE, &PLOT_COAST_CURVATURE_TITLE))
739 return false;
740 }
741
742 return true;
743 #else // #if defined(_SAGA_MSW) || defined(_SAGA_LINUX)
744 if( !bWriteVectorGIS(PLOT_COAST , (*m_pParameters)("COAST" )->asShapes()) )
745 return false;
746
747 if( !bWriteVectorGIS(PLOT_NORMALS , (*m_pParameters)("NORMALS" )->asShapes()) )
748 return false;
749
750 if( !bWriteVectorGIS(PLOT_CLIFF_TOP , (*m_pParameters)("CLIFF_TOP" )->asShapes()) )
751 return false;
752
753 if( !bWriteVectorGIS(PLOT_CLIFF_TOE , (*m_pParameters)("CLIFF_TOE" )->asShapes()) )
754 return false;
755
756 if( !bWriteVectorGIS(PLOT_COAST_POINT, (*m_pParameters)("COAST_POINT")->asShapes()) )
757 return false;
758
759 // These are optional (see delineation constructor)
760 if( (*m_pParameters)("INVALID_NORMALS")->asShapes() )
761 {
762 if( !bWriteVectorGIS(PLOT_INVALID_NORMALS, (*m_pParameters)("INVALID_NORMALS")->asShapes()) )
763 return false;
764 }
765
766 if( (*m_pParameters)("COAST_CURVATURE")->asShapes() )
767 {
768 if( !bWriteVectorGIS(PLOT_COAST_CURVATURE, (*m_pParameters)("COAST_CURVATURE")->asShapes()) )
769 return false;
770 }
771
772 return( true );
773 #endif // #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
774 }
775
776
777 /*==============================================================================================================================
778
779 Finds the max and min values in order to scale raster output if we cannot write floats or 32-bit integers
780
781 ==============================================================================================================================*/
GetRasterOutputMinMax(int const nDataItem,double & dMin,double & dMax)782 void CDelineation::GetRasterOutputMinMax(int const nDataItem, double& dMin, double& dMax)
783 {
784 // If this is a binary mask layer, we already know the max and min values
785 if ((nDataItem == PLOT_RASTER_COAST) ||
786 (nDataItem == PLOT_RASTER_NORMAL))
787
788 {
789 dMin = 0;
790 dMax = 1;
791
792 return;
793 }
794
795 // Not a binary mask layer, so we must find the max and min values
796 dMin = DBL_MAX;
797 dMax = DBL_MIN;
798
799 double dTmp;
800 for (int nY = 0; nY < m_nYGridMax; nY++)
801 {
802 for (int nX = 0; nX < m_nXGridMax; nX++)
803 {
804 switch (nDataItem)
805 {
806 case (PLOT_SEDIMENT_TOP_ELEV):
807 {
808 dTmp = m_pRasterGrid->pGetCell(nX, nY)->dGetSedimentTopElev();
809 break;
810 }
811 }
812
813 if (dTmp != DBL_NODATA)
814 {
815 if (dTmp > dMax)
816 dMax = dTmp;
817
818 if (dTmp < dMin)
819 dMin = dTmp;
820 }
821 }
822 }
823 }
824
825
826 /*==============================================================================================================================
827
828 Sets per-driver defaults for raster files created using GDAL
829
830 ===============================================================================================================================*/
SetRasterFileCreationDefaults(void)831 void CDelineation::SetRasterFileCreationDefaults(void)
832 {
833 #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
834 string
835 strDriver = strToLower(&m_strRasterGISOutFormat),
836 strComment = "Created by " + PROGNAME + " for " + PLATFORM + " " + strGetBuild() + " running on " + strGetComputerName();
837
838 // TODO Do these for all commonly-used file types
839 if (strDriver == "aaigrid")
840 {
841
842 }
843
844 else if (strDriver == "bmp")
845 {
846
847 }
848
849 else if (strDriver == "gtiff")
850 {
851 if (m_bWorldFile)
852 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "TFW", "YES");
853
854 // if (m_bCompressGTIFF)
855 // {
856 // m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "NUM_THREADS", "ALL_CPUS");
857 // m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "COMPRESS", "LZW");
858 // }
859 }
860
861 else if (strDriver == "hfa")
862 {
863 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "NBITS", "4");
864 }
865
866 else if (strDriver == "jpeg")
867 {
868 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "COMMENT", strComment.c_str());
869 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "QUALITY", "95");
870 }
871
872 else if (strDriver == "png")
873 {
874 if (m_bWorldFile)
875 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "WORLDFILE", "YES");
876
877 // m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "TITLE", "This is the title");
878 // m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "DESCRIPTION", "This is a description");
879 // m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "COPYRIGHT", "This is some copyright statement");
880 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "COMMENT", strComment.c_str());
881 m_papszGDALRasterOptions = CSLSetNameValue(m_papszGDALRasterOptions, "NBITS", "4");
882 }
883
884 else if (strDriver == "rst")
885 {
886
887 }
888
889 // for (int i = 0; m_papszGDALRasterOptions[i] != NULL; i++)
890 // cout << m_papszGDALRasterOptions[i] << endl;
891 // cout << endl;
892 #endif // #if !defined(_SAGA_MSW) && !defined(_SAGA_LINUX)
893 }
894
895
896
897
898