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