Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Natural Resources Canada's Geoid BYN file format
4 : * Purpose: Implementation of BYN format
5 : * Author: Ivan Lucena, ivan.lucena@outlook.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2018, Ivan Lucena
9 : * Copyright (c) 2018, Even Rouault
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_port.h"
31 : #include "byndataset.h"
32 : #include "rawdataset.h"
33 :
34 : #include "cpl_string.h"
35 : #include "gdal_frmts.h"
36 : #include "ogr_spatialref.h"
37 : #include "ogr_srs_api.h"
38 :
39 : #include <algorithm>
40 : #include <cstdlib>
41 : #include <limits>
42 :
43 : // Specification at
44 : // https://www.nrcan.gc.ca/sites/www.nrcan.gc.ca/files/earthsciences/pdf/gpshgrid_e.pdf
45 :
46 : const static BYNEllipsoids EllipsoidTable[] = {
47 : {"GRS80", 6378137.0, 298.257222101},
48 : {"WGS84", 6378137.0, 298.257223564},
49 : {"ALT1", 6378136.3, 298.256415099},
50 : {"GRS67", 6378160.0, 298.247167427},
51 : {"ELLIP1", 6378136.46, 298.256415099},
52 : {"ALT2", 6378136.3, 298.257},
53 : {"ELLIP2", 6378136.0, 298.257},
54 : {"CLARKE 1866", 6378206.4, 294.9786982}};
55 :
56 : /************************************************************************/
57 : /* BYNRasterBand() */
58 : /************************************************************************/
59 :
60 8 : BYNRasterBand::BYNRasterBand(GDALDataset *poDSIn, int nBandIn,
61 : VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
62 : int nPixelOffsetIn, int nLineOffsetIn,
63 8 : GDALDataType eDataTypeIn, int bNativeOrderIn)
64 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
65 : nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
66 8 : RawRasterBand::OwnFP::NO)
67 : {
68 8 : }
69 :
70 : /************************************************************************/
71 : /* ~BYNRasterBand() */
72 : /************************************************************************/
73 :
74 16 : BYNRasterBand::~BYNRasterBand()
75 : {
76 16 : }
77 :
78 : /************************************************************************/
79 : /* GetNoDataValue() */
80 : /************************************************************************/
81 :
82 6 : double BYNRasterBand::GetNoDataValue(int *pbSuccess)
83 : {
84 6 : if (pbSuccess)
85 5 : *pbSuccess = TRUE;
86 6 : int bSuccess = FALSE;
87 6 : double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
88 6 : if (bSuccess)
89 : {
90 2 : return dfNoData;
91 : }
92 4 : const double dfFactor =
93 4 : reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
94 4 : return eDataType == GDT_Int16 ? 32767.0 : 9999.0 * dfFactor;
95 : }
96 :
97 : /************************************************************************/
98 : /* GetScale() */
99 : /************************************************************************/
100 :
101 1 : double BYNRasterBand::GetScale(int *pbSuccess)
102 : {
103 1 : if (pbSuccess != nullptr)
104 1 : *pbSuccess = TRUE;
105 1 : const double dfFactor =
106 1 : reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
107 1 : return (dfFactor != 0.0) ? 1.0 / dfFactor : 0.0;
108 : }
109 :
110 : /************************************************************************/
111 : /* SetScale() */
112 : /************************************************************************/
113 :
114 1 : CPLErr BYNRasterBand::SetScale(double dfNewValue)
115 : {
116 1 : BYNDataset *poIDS = reinterpret_cast<BYNDataset *>(poDS);
117 1 : poIDS->hHeader.dfFactor = 1.0 / dfNewValue;
118 1 : return CE_None;
119 : }
120 :
121 : /************************************************************************/
122 : /* ==================================================================== */
123 : /* BYNDataset */
124 : /* ==================================================================== */
125 : /************************************************************************/
126 :
127 8 : BYNDataset::BYNDataset()
128 : : fpImage(nullptr), hHeader{0, 0, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0,
129 8 : 0, 0, 0, 0, 0, 0.0, 0.0, 0, 0, 0.0, 0}
130 : {
131 8 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
132 8 : adfGeoTransform[0] = 0.0;
133 8 : adfGeoTransform[1] = 1.0;
134 8 : adfGeoTransform[2] = 0.0;
135 8 : adfGeoTransform[3] = 0.0;
136 8 : adfGeoTransform[4] = 0.0;
137 8 : adfGeoTransform[5] = 1.0;
138 8 : }
139 :
140 : /************************************************************************/
141 : /* ~BYNDataset() */
142 : /************************************************************************/
143 :
144 16 : BYNDataset::~BYNDataset()
145 :
146 : {
147 8 : BYNDataset::Close();
148 16 : }
149 :
150 : /************************************************************************/
151 : /* Close() */
152 : /************************************************************************/
153 :
154 16 : CPLErr BYNDataset::Close()
155 : {
156 16 : CPLErr eErr = CE_None;
157 16 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
158 : {
159 8 : if (BYNDataset::FlushCache(true) != CE_None)
160 0 : eErr = CE_Failure;
161 :
162 8 : if (GetAccess() == GA_Update)
163 1 : UpdateHeader();
164 :
165 8 : if (fpImage != nullptr)
166 : {
167 8 : if (VSIFCloseL(fpImage) != 0)
168 : {
169 0 : eErr = CE_Failure;
170 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
171 : }
172 : }
173 :
174 8 : if (GDALPamDataset::Close() != CE_None)
175 0 : eErr = CE_Failure;
176 : }
177 16 : return eErr;
178 : }
179 :
180 : /************************************************************************/
181 : /* Identify() */
182 : /************************************************************************/
183 :
184 48609 : int BYNDataset::Identify(GDALOpenInfo *poOpenInfo)
185 :
186 : {
187 48609 : if (poOpenInfo->nHeaderBytes < BYN_HDR_SZ)
188 45554 : return FALSE;
189 :
190 : /* -------------------------------------------------------------------- */
191 : /* Check file extension (.byn/.err) */
192 : /* -------------------------------------------------------------------- */
193 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
194 3055 : const char *pszFileExtension = CPLGetExtension(poOpenInfo->pszFilename);
195 :
196 3056 : if (!EQUAL(pszFileExtension, "byn") && !EQUAL(pszFileExtension, "err"))
197 : {
198 3040 : return FALSE;
199 : }
200 : #endif
201 :
202 : /* -------------------------------------------------------------------- */
203 : /* Check some value's ranges on header */
204 : /* -------------------------------------------------------------------- */
205 :
206 16 : BYNHeader hHeader = {0, 0, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0,
207 : 0, 0, 0, 0, 0, 0.0, 0.0, 0, 0, 0.0, 0};
208 :
209 16 : buffer2header(poOpenInfo->pabyHeader, &hHeader);
210 :
211 16 : if (hHeader.nGlobal < 0 || hHeader.nGlobal > 1 || hHeader.nType < 0 ||
212 16 : hHeader.nType > 9 || (hHeader.nSizeOf != 2 && hHeader.nSizeOf != 4) ||
213 16 : hHeader.nVDatum < 0 || hHeader.nVDatum > 3 || hHeader.nDescrip < 0 ||
214 16 : hHeader.nDescrip > 3 || hHeader.nSubType < 0 || hHeader.nSubType > 9 ||
215 16 : hHeader.nDatum < 0 || hHeader.nDatum > 1 || hHeader.nEllipsoid < 0 ||
216 16 : hHeader.nEllipsoid > 7 || hHeader.nByteOrder < 0 ||
217 16 : hHeader.nByteOrder > 1 || hHeader.nScale < 0 || hHeader.nScale > 1)
218 0 : return FALSE;
219 :
220 : #if 0
221 : // We have disabled those checks as invalid values are often found in some
222 : // datasets, such as http://s3.microsurvey.com/os/fieldgenius/geoids/Lithuania.zip
223 : // We don't use those fields, so we may just ignore them.
224 : if((hHeader.nTideSys < 0 || hHeader.nTideSys > 2 ||
225 : hHeader.nPtType < 0 || hHeader.nPtType > 1 ))
226 : {
227 : // Some datasets use 0xCC as a marker for invalidity for
228 : // records starting from Geopotential Wo
229 : for( int i = 52; i < 78; i++ )
230 : {
231 : if( poOpenInfo->pabyHeader[i] != 0xCC )
232 : return FALSE;
233 : }
234 : }
235 : #endif
236 :
237 16 : if (hHeader.nScale == 0)
238 : {
239 32 : if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
240 32 : (hHeader.nDLat / 2)) > BYN_MAX_LAT) ||
241 16 : (std::abs(static_cast<GIntBig>(hHeader.nNorth) +
242 32 : (hHeader.nDLat / 2)) > BYN_MAX_LAT) ||
243 16 : (std::abs(static_cast<GIntBig>(hHeader.nWest) -
244 48 : (hHeader.nDLon / 2)) > BYN_MAX_LON) ||
245 16 : (std::abs(static_cast<GIntBig>(hHeader.nEast) +
246 16 : (hHeader.nDLon / 2)) > BYN_MAX_LON))
247 0 : return FALSE;
248 : }
249 : else
250 : {
251 0 : if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
252 0 : (hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) ||
253 0 : (std::abs(static_cast<GIntBig>(hHeader.nNorth) +
254 0 : (hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) ||
255 0 : (std::abs(static_cast<GIntBig>(hHeader.nWest) -
256 0 : (hHeader.nDLon / 2)) > BYN_MAX_LON_SCL) ||
257 0 : (std::abs(static_cast<GIntBig>(hHeader.nEast) +
258 0 : (hHeader.nDLon / 2)) > BYN_MAX_LON_SCL))
259 0 : return FALSE;
260 : }
261 :
262 16 : return TRUE;
263 : }
264 :
265 : /************************************************************************/
266 : /* Open() */
267 : /************************************************************************/
268 :
269 8 : GDALDataset *BYNDataset::Open(GDALOpenInfo *poOpenInfo)
270 :
271 : {
272 8 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
273 0 : return nullptr;
274 :
275 : /* -------------------------------------------------------------------- */
276 : /* Create a corresponding GDALDataset. */
277 : /* -------------------------------------------------------------------- */
278 :
279 16 : auto poDS = std::make_unique<BYNDataset>();
280 :
281 8 : poDS->eAccess = poOpenInfo->eAccess;
282 8 : std::swap(poDS->fpImage, poOpenInfo->fpL);
283 :
284 : /* -------------------------------------------------------------------- */
285 : /* Read the header. */
286 : /* -------------------------------------------------------------------- */
287 :
288 8 : buffer2header(poOpenInfo->pabyHeader, &poDS->hHeader);
289 :
290 : /********************************/
291 : /* Scale boundaries and spacing */
292 : /********************************/
293 :
294 8 : double dfSouth = poDS->hHeader.nSouth;
295 8 : double dfNorth = poDS->hHeader.nNorth;
296 8 : double dfWest = poDS->hHeader.nWest;
297 8 : double dfEast = poDS->hHeader.nEast;
298 8 : double dfDLat = poDS->hHeader.nDLat;
299 8 : double dfDLon = poDS->hHeader.nDLon;
300 :
301 8 : if (poDS->hHeader.nScale == 1)
302 : {
303 0 : dfSouth *= BYN_SCALE;
304 0 : dfNorth *= BYN_SCALE;
305 0 : dfWest *= BYN_SCALE;
306 0 : dfEast *= BYN_SCALE;
307 0 : dfDLat *= BYN_SCALE;
308 0 : dfDLon *= BYN_SCALE;
309 : }
310 :
311 : /******************************/
312 : /* Calculate rows and columns */
313 : /******************************/
314 :
315 8 : double dfXSize = -1;
316 8 : double dfYSize = -1;
317 :
318 8 : poDS->nRasterXSize = -1;
319 8 : poDS->nRasterYSize = -1;
320 :
321 8 : if (dfDLat != 0.0 && dfDLon != 0.0)
322 : {
323 8 : dfXSize = ((dfEast - dfWest + 1.0) / dfDLon) + 1.0;
324 8 : dfYSize = ((dfNorth - dfSouth + 1.0) / dfDLat) + 1.0;
325 : }
326 :
327 8 : if (dfXSize > 0.0 && dfXSize < std::numeric_limits<double>::max() &&
328 16 : dfYSize > 0.0 && dfYSize < std::numeric_limits<double>::max())
329 : {
330 8 : poDS->nRasterXSize = static_cast<GInt32>(dfXSize);
331 8 : poDS->nRasterYSize = static_cast<GInt32>(dfYSize);
332 : }
333 :
334 8 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
335 : {
336 0 : return nullptr;
337 : }
338 :
339 : /*****************************/
340 : /* Build GeoTransform matrix */
341 : /*****************************/
342 :
343 8 : poDS->adfGeoTransform[0] = (dfWest - (dfDLon / 2.0)) / 3600.0;
344 8 : poDS->adfGeoTransform[1] = dfDLon / 3600.0;
345 8 : poDS->adfGeoTransform[2] = 0.0;
346 8 : poDS->adfGeoTransform[3] = (dfNorth + (dfDLat / 2.0)) / 3600.0;
347 8 : poDS->adfGeoTransform[4] = 0.0;
348 8 : poDS->adfGeoTransform[5] = -1 * dfDLat / 3600.0;
349 :
350 : /*********************/
351 : /* Set data type */
352 : /*********************/
353 :
354 8 : GDALDataType eDT = GDT_Unknown;
355 :
356 8 : if (poDS->hHeader.nSizeOf == 2)
357 0 : eDT = GDT_Int16;
358 8 : else if (poDS->hHeader.nSizeOf == 4)
359 8 : eDT = GDT_Int32;
360 : else
361 : {
362 0 : return nullptr;
363 : }
364 :
365 : /* -------------------------------------------------------------------- */
366 : /* Create band information object. */
367 : /* -------------------------------------------------------------------- */
368 :
369 8 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
370 :
371 8 : int bIsLSB = poDS->hHeader.nByteOrder == 1 ? 1 : 0;
372 :
373 : auto poBand = std::make_unique<BYNRasterBand>(
374 8 : poDS.get(), 1, poDS->fpImage, BYN_HDR_SZ, nDTSize,
375 24 : poDS->nRasterXSize * nDTSize, eDT, CPL_IS_LSB == bIsLSB);
376 8 : if (!poBand->IsValid())
377 0 : return nullptr;
378 8 : poDS->SetBand(1, std::move(poBand));
379 :
380 : /* -------------------------------------------------------------------- */
381 : /* Initialize any PAM information. */
382 : /* -------------------------------------------------------------------- */
383 :
384 8 : poDS->SetDescription(poOpenInfo->pszFilename);
385 8 : poDS->TryLoadXML();
386 :
387 : /* -------------------------------------------------------------------- */
388 : /* Check for overviews. */
389 : /* -------------------------------------------------------------------- */
390 :
391 8 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
392 :
393 8 : return poDS.release();
394 : }
395 :
396 : /************************************************************************/
397 : /* GetGeoTransform() */
398 : /************************************************************************/
399 :
400 2 : CPLErr BYNDataset::GetGeoTransform(double *padfTransform)
401 :
402 : {
403 2 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
404 2 : return CE_None;
405 : }
406 :
407 : /************************************************************************/
408 : /* SetGeoTransform() */
409 : /************************************************************************/
410 :
411 1 : CPLErr BYNDataset::SetGeoTransform(double *padfTransform)
412 :
413 : {
414 1 : if (padfTransform[2] != 0.0 || padfTransform[4] != 0.0)
415 : {
416 0 : CPLError(CE_Failure, CPLE_AppDefined,
417 : "Attempt to write skewed or rotated geotransform to byn.");
418 0 : return CE_Failure;
419 : }
420 1 : memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
421 1 : return CE_None;
422 : }
423 :
424 : /************************************************************************/
425 : /* GetSpatialRef() */
426 : /************************************************************************/
427 :
428 1 : const OGRSpatialReference *BYNDataset::GetSpatialRef() const
429 :
430 : {
431 1 : if (!m_oSRS.IsEmpty())
432 0 : return &m_oSRS;
433 :
434 : /* Try to use a prefefined EPSG compound CS */
435 :
436 1 : if (hHeader.nDatum == 1 && hHeader.nVDatum == 2)
437 : {
438 0 : m_oSRS.importFromEPSG(BYN_DATUM_1_VDATUM_2);
439 0 : return &m_oSRS;
440 : }
441 :
442 : /* Build the GEOGCS based on Datum ( or Ellipsoid )*/
443 :
444 1 : bool bNoGeogCS = false;
445 :
446 1 : if (hHeader.nDatum == 0)
447 1 : m_oSRS.importFromEPSG(BYN_DATUM_0);
448 0 : else if (hHeader.nDatum == 1)
449 0 : m_oSRS.importFromEPSG(BYN_DATUM_1);
450 : else
451 : {
452 : /* Build GEOGCS based on Ellipsoid (Table 3) */
453 :
454 0 : if (hHeader.nEllipsoid > -1 &&
455 0 : hHeader.nEllipsoid <
456 : static_cast<GInt16>(CPL_ARRAYSIZE(EllipsoidTable)))
457 0 : m_oSRS.SetGeogCS(
458 0 : CPLSPrintf("BYN Ellipsoid(%d)", hHeader.nEllipsoid),
459 0 : "Unspecified", EllipsoidTable[hHeader.nEllipsoid].pszName,
460 0 : EllipsoidTable[hHeader.nEllipsoid].dfSemiMajor,
461 0 : EllipsoidTable[hHeader.nEllipsoid].dfInvFlattening);
462 : else
463 0 : bNoGeogCS = true;
464 : }
465 :
466 : /* Build the VERT_CS based on VDatum */
467 :
468 2 : OGRSpatialReference oSRSComp;
469 2 : OGRSpatialReference oSRSVert;
470 :
471 1 : int nVertCS = 0;
472 :
473 1 : if (hHeader.nVDatum == 1)
474 0 : nVertCS = BYN_VDATUM_1;
475 1 : else if (hHeader.nVDatum == 2)
476 1 : nVertCS = BYN_VDATUM_2;
477 0 : else if (hHeader.nVDatum == 3)
478 0 : nVertCS = BYN_VDATUM_3;
479 : else
480 : {
481 : /* Return GEOGCS ( .err files ) */
482 :
483 0 : if (bNoGeogCS)
484 0 : return nullptr;
485 :
486 0 : return &m_oSRS;
487 : }
488 :
489 1 : oSRSVert.importFromEPSG(nVertCS);
490 :
491 : /* Create CPMPD_CS with GEOGCS and VERT_CS */
492 :
493 1 : if (oSRSComp.SetCompoundCS(CPLSPrintf("BYN Datum(%d) & VDatum(%d)",
494 1 : hHeader.nDatum, hHeader.nDatum),
495 2 : &m_oSRS, &oSRSVert) == CE_None)
496 : {
497 : /* Return COMPD_CS with GEOGCS and VERT_CS */
498 :
499 1 : m_oSRS = std::move(oSRSComp);
500 1 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
501 1 : return &m_oSRS;
502 : }
503 :
504 0 : return nullptr;
505 : }
506 :
507 : /************************************************************************/
508 : /* SetSpatialRef() */
509 : /************************************************************************/
510 :
511 1 : CPLErr BYNDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
512 :
513 : {
514 1 : if (poSRS == nullptr)
515 0 : return CE_Failure;
516 :
517 1 : m_oSRS = *poSRS;
518 :
519 : /* Try to recognize prefefined EPSG compound CS */
520 :
521 1 : if (poSRS->IsCompound())
522 : {
523 1 : const char *pszAuthName = poSRS->GetAuthorityName("COMPD_CS");
524 1 : const char *pszAuthCode = poSRS->GetAuthorityCode("COMPD_CS");
525 :
526 1 : if (pszAuthName != nullptr && pszAuthCode != nullptr &&
527 0 : EQUAL(pszAuthName, "EPSG") &&
528 0 : atoi(pszAuthCode) == BYN_DATUM_1_VDATUM_2)
529 : {
530 0 : hHeader.nVDatum = 2;
531 0 : hHeader.nDatum = 1;
532 0 : return CE_None;
533 : }
534 : }
535 :
536 1 : OGRSpatialReference oSRSTemp;
537 :
538 : /* Try to match GEOGCS */
539 :
540 1 : if (poSRS->IsGeographic())
541 : {
542 1 : oSRSTemp.importFromEPSG(BYN_DATUM_0);
543 1 : if (poSRS->IsSameGeogCS(&oSRSTemp))
544 1 : hHeader.nDatum = 0;
545 : else
546 : {
547 0 : oSRSTemp.importFromEPSG(BYN_DATUM_1);
548 0 : if (poSRS->IsSameGeogCS(&oSRSTemp))
549 0 : hHeader.nDatum = 1;
550 : }
551 : }
552 :
553 : /* Try to match VERT_CS */
554 :
555 1 : if (poSRS->IsVertical())
556 : {
557 1 : oSRSTemp.importFromEPSG(BYN_VDATUM_1);
558 1 : if (poSRS->IsSameVertCS(&oSRSTemp))
559 0 : hHeader.nVDatum = 1;
560 : else
561 : {
562 1 : oSRSTemp.importFromEPSG(BYN_VDATUM_2);
563 1 : if (poSRS->IsSameVertCS(&oSRSTemp))
564 1 : hHeader.nVDatum = 2;
565 : else
566 : {
567 0 : oSRSTemp.importFromEPSG(BYN_VDATUM_3);
568 0 : if (poSRS->IsSameVertCS(&oSRSTemp))
569 0 : hHeader.nVDatum = 3;
570 : }
571 : }
572 : }
573 :
574 1 : return CE_None;
575 : }
576 :
577 : /*----------------------------------------------------------------------*/
578 : /* buffer2header() */
579 : /*----------------------------------------------------------------------*/
580 :
581 24 : void BYNDataset::buffer2header(const GByte *pabyBuf, BYNHeader *pohHeader)
582 :
583 : {
584 24 : memcpy(&pohHeader->nSouth, pabyBuf, 4);
585 24 : memcpy(&pohHeader->nNorth, pabyBuf + 4, 4);
586 24 : memcpy(&pohHeader->nWest, pabyBuf + 8, 4);
587 24 : memcpy(&pohHeader->nEast, pabyBuf + 12, 4);
588 24 : memcpy(&pohHeader->nDLat, pabyBuf + 16, 2);
589 24 : memcpy(&pohHeader->nDLon, pabyBuf + 18, 2);
590 24 : memcpy(&pohHeader->nGlobal, pabyBuf + 20, 2);
591 24 : memcpy(&pohHeader->nType, pabyBuf + 22, 2);
592 24 : memcpy(&pohHeader->dfFactor, pabyBuf + 24, 8);
593 24 : memcpy(&pohHeader->nSizeOf, pabyBuf + 32, 2);
594 24 : memcpy(&pohHeader->nVDatum, pabyBuf + 34, 2);
595 24 : memcpy(&pohHeader->nDescrip, pabyBuf + 40, 2);
596 24 : memcpy(&pohHeader->nSubType, pabyBuf + 42, 2);
597 24 : memcpy(&pohHeader->nDatum, pabyBuf + 44, 2);
598 24 : memcpy(&pohHeader->nEllipsoid, pabyBuf + 46, 2);
599 24 : memcpy(&pohHeader->nByteOrder, pabyBuf + 48, 2);
600 24 : memcpy(&pohHeader->nScale, pabyBuf + 50, 2);
601 24 : memcpy(&pohHeader->dfWo, pabyBuf + 52, 8);
602 24 : memcpy(&pohHeader->dfGM, pabyBuf + 60, 8);
603 24 : memcpy(&pohHeader->nTideSys, pabyBuf + 68, 2);
604 24 : memcpy(&pohHeader->nRealiz, pabyBuf + 70, 2);
605 24 : memcpy(&pohHeader->dEpoch, pabyBuf + 72, 4);
606 24 : memcpy(&pohHeader->nPtType, pabyBuf + 76, 2);
607 :
608 : #if defined(CPL_MSB)
609 : CPL_LSBPTR32(&pohHeader->nSouth);
610 : CPL_LSBPTR32(&pohHeader->nNorth);
611 : CPL_LSBPTR32(&pohHeader->nWest);
612 : CPL_LSBPTR32(&pohHeader->nEast);
613 : CPL_LSBPTR16(&pohHeader->nDLat);
614 : CPL_LSBPTR16(&pohHeader->nDLon);
615 : CPL_LSBPTR16(&pohHeader->nGlobal);
616 : CPL_LSBPTR16(&pohHeader->nType);
617 : CPL_LSBPTR64(&pohHeader->dfFactor);
618 : CPL_LSBPTR16(&pohHeader->nSizeOf);
619 : CPL_LSBPTR16(&pohHeader->nVDatum);
620 : CPL_LSBPTR16(&pohHeader->nDescrip);
621 : CPL_LSBPTR16(&pohHeader->nSubType);
622 : CPL_LSBPTR16(&pohHeader->nDatum);
623 : CPL_LSBPTR16(&pohHeader->nEllipsoid);
624 : CPL_LSBPTR16(&pohHeader->nByteOrder);
625 : CPL_LSBPTR16(&pohHeader->nScale);
626 : CPL_LSBPTR64(&pohHeader->dfWo);
627 : CPL_LSBPTR64(&pohHeader->dfGM);
628 : CPL_LSBPTR16(&pohHeader->nTideSys);
629 : CPL_LSBPTR16(&pohHeader->nRealiz);
630 : CPL_LSBPTR32(&pohHeader->dEpoch);
631 : CPL_LSBPTR16(&pohHeader->nPtType);
632 : #endif
633 :
634 : #if DEBUG
635 24 : CPLDebug("BYN", "South = %d", pohHeader->nSouth);
636 24 : CPLDebug("BYN", "North = %d", pohHeader->nNorth);
637 24 : CPLDebug("BYN", "West = %d", pohHeader->nWest);
638 24 : CPLDebug("BYN", "East = %d", pohHeader->nEast);
639 24 : CPLDebug("BYN", "DLat = %d", pohHeader->nDLat);
640 24 : CPLDebug("BYN", "DLon = %d", pohHeader->nDLon);
641 24 : CPLDebug("BYN", "DGlobal = %d", pohHeader->nGlobal);
642 24 : CPLDebug("BYN", "DType = %d", pohHeader->nType);
643 24 : CPLDebug("BYN", "Factor = %f", pohHeader->dfFactor);
644 24 : CPLDebug("BYN", "SizeOf = %d", pohHeader->nSizeOf);
645 24 : CPLDebug("BYN", "VDatum = %d", pohHeader->nVDatum);
646 24 : CPLDebug("BYN", "Data = %d", pohHeader->nDescrip);
647 24 : CPLDebug("BYN", "SubType = %d", pohHeader->nSubType);
648 24 : CPLDebug("BYN", "Datum = %d", pohHeader->nDatum);
649 24 : CPLDebug("BYN", "Ellipsoid = %d", pohHeader->nEllipsoid);
650 24 : CPLDebug("BYN", "ByteOrder = %d", pohHeader->nByteOrder);
651 24 : CPLDebug("BYN", "Scale = %d", pohHeader->nScale);
652 24 : CPLDebug("BYN", "Wo = %f", pohHeader->dfWo);
653 24 : CPLDebug("BYN", "GM = %f", pohHeader->dfGM);
654 24 : CPLDebug("BYN", "TideSystem = %d", pohHeader->nTideSys);
655 24 : CPLDebug("BYN", "RefRealzation = %d", pohHeader->nRealiz);
656 24 : CPLDebug("BYN", "Epoch = %f", pohHeader->dEpoch);
657 24 : CPLDebug("BYN", "PtType = %d", pohHeader->nPtType);
658 : #endif
659 24 : }
660 :
661 : /*----------------------------------------------------------------------*/
662 : /* header2buffer() */
663 : /*----------------------------------------------------------------------*/
664 :
665 2 : void BYNDataset::header2buffer(const BYNHeader *pohHeader, GByte *pabyBuf)
666 :
667 : {
668 2 : memcpy(pabyBuf, &pohHeader->nSouth, 4);
669 2 : memcpy(pabyBuf + 4, &pohHeader->nNorth, 4);
670 2 : memcpy(pabyBuf + 8, &pohHeader->nWest, 4);
671 2 : memcpy(pabyBuf + 12, &pohHeader->nEast, 4);
672 2 : memcpy(pabyBuf + 16, &pohHeader->nDLat, 2);
673 2 : memcpy(pabyBuf + 18, &pohHeader->nDLon, 2);
674 2 : memcpy(pabyBuf + 20, &pohHeader->nGlobal, 2);
675 2 : memcpy(pabyBuf + 22, &pohHeader->nType, 2);
676 2 : memcpy(pabyBuf + 24, &pohHeader->dfFactor, 8);
677 2 : memcpy(pabyBuf + 32, &pohHeader->nSizeOf, 2);
678 2 : memcpy(pabyBuf + 34, &pohHeader->nVDatum, 2);
679 2 : memcpy(pabyBuf + 40, &pohHeader->nDescrip, 2);
680 2 : memcpy(pabyBuf + 42, &pohHeader->nSubType, 2);
681 2 : memcpy(pabyBuf + 44, &pohHeader->nDatum, 2);
682 2 : memcpy(pabyBuf + 46, &pohHeader->nEllipsoid, 2);
683 2 : memcpy(pabyBuf + 48, &pohHeader->nByteOrder, 2);
684 2 : memcpy(pabyBuf + 50, &pohHeader->nScale, 2);
685 2 : memcpy(pabyBuf + 52, &pohHeader->dfWo, 8);
686 2 : memcpy(pabyBuf + 60, &pohHeader->dfGM, 8);
687 2 : memcpy(pabyBuf + 68, &pohHeader->nTideSys, 2);
688 2 : memcpy(pabyBuf + 70, &pohHeader->nRealiz, 2);
689 2 : memcpy(pabyBuf + 72, &pohHeader->dEpoch, 4);
690 2 : memcpy(pabyBuf + 76, &pohHeader->nPtType, 2);
691 :
692 : #if defined(CPL_MSB)
693 : CPL_LSBPTR32(pabyBuf);
694 : CPL_LSBPTR32(pabyBuf + 4);
695 : CPL_LSBPTR32(pabyBuf + 8);
696 : CPL_LSBPTR32(pabyBuf + 12);
697 : CPL_LSBPTR16(pabyBuf + 16);
698 : CPL_LSBPTR16(pabyBuf + 18);
699 : CPL_LSBPTR16(pabyBuf + 20);
700 : CPL_LSBPTR16(pabyBuf + 22);
701 : CPL_LSBPTR64(pabyBuf + 24);
702 : CPL_LSBPTR16(pabyBuf + 32);
703 : CPL_LSBPTR16(pabyBuf + 34);
704 : CPL_LSBPTR16(pabyBuf + 40);
705 : CPL_LSBPTR16(pabyBuf + 42);
706 : CPL_LSBPTR16(pabyBuf + 44);
707 : CPL_LSBPTR16(pabyBuf + 46);
708 : CPL_LSBPTR16(pabyBuf + 48);
709 : CPL_LSBPTR16(pabyBuf + 50);
710 : CPL_LSBPTR64(pabyBuf + 52);
711 : CPL_LSBPTR64(pabyBuf + 60);
712 : CPL_LSBPTR16(pabyBuf + 68);
713 : CPL_LSBPTR16(pabyBuf + 70);
714 : CPL_LSBPTR32(pabyBuf + 72);
715 : CPL_LSBPTR16(pabyBuf + 76);
716 : #endif
717 2 : }
718 :
719 : /************************************************************************/
720 : /* Create() */
721 : /************************************************************************/
722 :
723 50 : GDALDataset *BYNDataset::Create(const char *pszFilename, int nXSize, int nYSize,
724 : int /* nBands */, GDALDataType eType,
725 : char ** /* papszOptions */)
726 : {
727 50 : if (eType != GDT_Int16 && eType != GDT_Int32)
728 : {
729 43 : CPLError(CE_Failure, CPLE_AppDefined,
730 : "Attempt to create byn file with unsupported data type '%s'.",
731 : GDALGetDataTypeName(eType));
732 43 : return nullptr;
733 : }
734 :
735 : /* -------------------------------------------------------------------- */
736 : /* Check file extension (.byn/.err) */
737 : /* -------------------------------------------------------------------- */
738 :
739 7 : char *pszFileExtension = CPLStrdup(CPLGetExtension(pszFilename));
740 :
741 7 : if (!EQUAL(pszFileExtension, "byn") && !EQUAL(pszFileExtension, "err"))
742 : {
743 6 : CPLError(
744 : CE_Failure, CPLE_AppDefined,
745 : "Attempt to create byn file with extension other than byn/err.");
746 6 : CPLFree(pszFileExtension);
747 6 : return nullptr;
748 : }
749 :
750 1 : CPLFree(pszFileExtension);
751 :
752 : /* -------------------------------------------------------------------- */
753 : /* Try to create the file. */
754 : /* -------------------------------------------------------------------- */
755 :
756 1 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb+");
757 1 : if (fp == nullptr)
758 : {
759 0 : CPLError(CE_Failure, CPLE_OpenFailed,
760 : "Attempt to create file `%s' failed.\n", pszFilename);
761 0 : return nullptr;
762 : }
763 :
764 : /* -------------------------------------------------------------------- */
765 : /* Write an empty header. */
766 : /* -------------------------------------------------------------------- */
767 :
768 1 : GByte abyBuf[BYN_HDR_SZ] = {'\0'};
769 :
770 : /* Load header with some commum values */
771 :
772 1 : BYNHeader hHeader = {0, 0, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0,
773 : 0, 0, 0, 0, 0, 0.0, 0.0, 0, 0, 0.0, 0};
774 :
775 : /* Set temporary values on header */
776 :
777 1 : hHeader.nSouth = 0;
778 1 : hHeader.nNorth = nYSize - 2;
779 1 : hHeader.nWest = 0;
780 1 : hHeader.nEast = nXSize - 2;
781 1 : hHeader.nDLat = 1;
782 1 : hHeader.nDLon = 1;
783 1 : hHeader.nSizeOf = static_cast<GInt16>(GDALGetDataTypeSizeBytes(eType));
784 :
785 : /* Prepare buffer for writing */
786 :
787 1 : header2buffer(&hHeader, abyBuf);
788 :
789 : /* Write initial header */
790 :
791 1 : VSIFWriteL(abyBuf, BYN_HDR_SZ, 1, fp);
792 1 : VSIFCloseL(fp);
793 :
794 1 : return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
795 : }
796 :
797 : /************************************************************************/
798 : /* UpdateHeader() */
799 : /************************************************************************/
800 :
801 1 : void BYNDataset::UpdateHeader()
802 : {
803 1 : double dfDLon = (adfGeoTransform[1] * 3600.0);
804 1 : double dfDLat = (adfGeoTransform[5] * 3600.0 * -1);
805 1 : double dfWest = (adfGeoTransform[0] * 3600.0) + (dfDLon / 2);
806 1 : double dfNorth = (adfGeoTransform[3] * 3600.0) - (dfDLat / 2);
807 1 : double dfSouth = dfNorth - ((nRasterYSize - 1) * dfDLat);
808 1 : double dfEast = dfWest + ((nRasterXSize - 1) * dfDLon);
809 :
810 1 : if (hHeader.nScale == 1)
811 : {
812 0 : dfSouth /= BYN_SCALE;
813 0 : dfNorth /= BYN_SCALE;
814 0 : dfWest /= BYN_SCALE;
815 0 : dfEast /= BYN_SCALE;
816 0 : dfDLat /= BYN_SCALE;
817 0 : dfDLon /= BYN_SCALE;
818 : }
819 :
820 1 : hHeader.nSouth = static_cast<GInt32>(dfSouth);
821 1 : hHeader.nNorth = static_cast<GInt32>(dfNorth);
822 1 : hHeader.nWest = static_cast<GInt32>(dfWest);
823 1 : hHeader.nEast = static_cast<GInt32>(dfEast);
824 1 : hHeader.nDLat = static_cast<GInt16>(dfDLat);
825 1 : hHeader.nDLon = static_cast<GInt16>(dfDLon);
826 :
827 : GByte abyBuf[BYN_HDR_SZ];
828 :
829 1 : header2buffer(&hHeader, abyBuf);
830 :
831 1 : const char *pszValue = GetMetadataItem("GLOBAL");
832 1 : if (pszValue != nullptr)
833 0 : hHeader.nGlobal = static_cast<GInt16>(atoi(pszValue));
834 :
835 1 : pszValue = GetMetadataItem("TYPE");
836 1 : if (pszValue != nullptr)
837 0 : hHeader.nType = static_cast<GInt16>(atoi(pszValue));
838 :
839 1 : pszValue = GetMetadataItem("DESCRIPTION");
840 1 : if (pszValue != nullptr)
841 0 : hHeader.nDescrip = static_cast<GInt16>(atoi(pszValue));
842 :
843 1 : pszValue = GetMetadataItem("SUBTYPE");
844 1 : if (pszValue != nullptr)
845 0 : hHeader.nSubType = static_cast<GInt16>(atoi(pszValue));
846 :
847 1 : pszValue = GetMetadataItem("WO");
848 1 : if (pszValue != nullptr)
849 0 : hHeader.dfWo = CPLAtof(pszValue);
850 :
851 1 : pszValue = GetMetadataItem("GM");
852 1 : if (pszValue != nullptr)
853 0 : hHeader.dfGM = CPLAtof(pszValue);
854 :
855 1 : pszValue = GetMetadataItem("TIDESYSTEM");
856 1 : if (pszValue != nullptr)
857 0 : hHeader.nTideSys = static_cast<GInt16>(atoi(pszValue));
858 :
859 1 : pszValue = GetMetadataItem("REALIZATION");
860 1 : if (pszValue != nullptr)
861 0 : hHeader.nRealiz = static_cast<GInt16>(atoi(pszValue));
862 :
863 1 : pszValue = GetMetadataItem("EPOCH");
864 1 : if (pszValue != nullptr)
865 0 : hHeader.dEpoch = static_cast<float>(CPLAtof(pszValue));
866 :
867 1 : pszValue = GetMetadataItem("PTTYPE");
868 1 : if (pszValue != nullptr)
869 0 : hHeader.nPtType = static_cast<GInt16>(atoi(pszValue));
870 :
871 1 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
872 1 : CPL_IGNORE_RET_VAL(VSIFWriteL(abyBuf, BYN_HDR_SZ, 1, fpImage));
873 :
874 : /* GDALPam metadata update */
875 :
876 1 : SetMetadataItem("GLOBAL", CPLSPrintf("%d", hHeader.nGlobal), "BYN");
877 1 : SetMetadataItem("TYPE", CPLSPrintf("%d", hHeader.nType), "BYN");
878 1 : SetMetadataItem("DESCRIPTION", CPLSPrintf("%d", hHeader.nDescrip), "BYN");
879 1 : SetMetadataItem("SUBTYPE", CPLSPrintf("%d", hHeader.nSubType), "BYN");
880 1 : SetMetadataItem("WO", CPLSPrintf("%g", hHeader.dfWo), "BYN");
881 1 : SetMetadataItem("GM", CPLSPrintf("%g", hHeader.dfGM), "BYN");
882 1 : SetMetadataItem("TIDESYSTEM", CPLSPrintf("%d", hHeader.nTideSys), "BYN");
883 1 : SetMetadataItem("REALIZATION", CPLSPrintf("%d", hHeader.nRealiz), "BYN");
884 1 : SetMetadataItem("EPOCH", CPLSPrintf("%g", hHeader.dEpoch), "BYN");
885 1 : SetMetadataItem("PTTYPE", CPLSPrintf("%d", hHeader.nPtType), "BYN");
886 1 : }
887 :
888 : /************************************************************************/
889 : /* GDALRegister_BYN() */
890 : /************************************************************************/
891 :
892 1511 : void GDALRegister_BYN()
893 :
894 : {
895 1511 : if (GDALGetDriverByName("BYN") != nullptr)
896 295 : return;
897 :
898 1216 : GDALDriver *poDriver = new GDALDriver();
899 :
900 1216 : poDriver->SetDescription("BYN");
901 1216 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
902 1216 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
903 1216 : "Natural Resources Canada's Geoid");
904 1216 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "byn err");
905 1216 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
906 1216 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/byn.html");
907 1216 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Int16 Int32");
908 :
909 1216 : poDriver->pfnOpen = BYNDataset::Open;
910 1216 : poDriver->pfnIdentify = BYNDataset::Identify;
911 1216 : poDriver->pfnCreate = BYNDataset::Create;
912 :
913 1216 : GetGDALDriverManager()->RegisterDriver(poDriver);
914 : }
|
