Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Implementation for NSIDC binary format.
5 : * Author: Michael Sumner, mdsumner@gmail.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2022, Michael Sumner
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_string.h"
14 : #include "gdal_frmts.h"
15 : #include "ogr_spatialref.h"
16 : #include "rawdataset.h"
17 :
18 : #include <cmath>
19 : #include <algorithm>
20 :
21 : /***********************************************************************/
22 : /* ====================================================================*/
23 : /* NSIDCbinDataset */
24 : /* ====================================================================*/
25 : /***********************************************************************/
26 :
27 : class NSIDCbinDataset final : public GDALPamDataset
28 : {
29 : friend class NSIDCbinRasterBand;
30 :
31 : struct NSIDCbinHeader
32 : {
33 :
34 : // page 7, User Guide https://nsidc.org/data/nsidc-0051
35 : // 1.3.2 File Contents
36 : // The file format consists of a 300-byte descriptive header followed by
37 : // a two-dimensional array of one-byte values containing the data. The
38 : // file header is composed of:
39 : // - a 21-element array of 6-byte character strings that contain
40 : // information such as polar stereographic grid characteristics
41 : // - a 24-byte character string containing the file name
42 : // - a 80-character string containing an optional image title
43 : // - a 70-byte character string containing ancillary information such as
44 : // data origin, data set creation date, etc. For compatibility with ANSI
45 : // C, IDL, and other languages, character strings are terminated with a
46 : // NULL byte.
47 : // Example file:
48 : // ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/south/daily/2010/nt_20100918_f17_v1.1_s.bin
49 :
50 : char missing_int[6] = {0}; // "00255"
51 : char columns[6] = {0}; // " 316"
52 : char rows[6] = {0}; // " 332"
53 : char internal1[6] = {0}; // "1.799"
54 : char latitude[6] = {0}; // "-51.3"
55 : char greenwich[6] = {0}; // "270.0"
56 : char internal2[6] = {0}; // "558.4"
57 : char jpole[6] = {0}; // "158.0"
58 : char ipole[6] = {0}; // "174.0"
59 : char instrument[6] = {0}; // "SSMIS"
60 : char data_descriptors[6] = {0}; // "17 cn"
61 : char julian_start[6] = {0}; // " 261"
62 : char hour_start[6] = {0}; // "-9999"
63 : char minute_start[6] = {0}; // "-9999"
64 : char julian_end[6] = {0}; // " 261"
65 : char hour_end[6] = {0}; // "-9999"
66 : char minute_end[6] = {0}; // "-9999"
67 : char year[6] = {0}; // " 2010"
68 : char julian[6] = {0}; // " 261"
69 : char channel[6] = {0}; // " 000"
70 : char scaling[6] = {0}; // "00250"
71 :
72 : // 121-126 Integer scaling factor
73 : // 127-150 24-character file name (without file-name extension)
74 : // 151-230 80-character image title
75 : // 231-300 70-character data information (creation date, data source,
76 : // etc.)
77 : char filename[24] = {0}; // " nt_20100918_f17_v1.1_s"
78 : char imagetitle[80] = {0}; // "ANTARCTIC SMMR TOTAL ICE CONCENTRATION
79 : // NIMBUSN07 DAY 299 10/26/1978"
80 : char data_information[70] = {0}; // "ANTARCTIC SMMR ONSSMIGRID CON
81 : // Coast253Pole251Land254 06/27/1996"
82 : };
83 :
84 : VSILFILE *fp = nullptr;
85 : CPLString osSRS{};
86 : NSIDCbinHeader sHeader{};
87 :
88 : double adfGeoTransform[6];
89 : CPL_DISALLOW_COPY_ASSIGN(NSIDCbinDataset)
90 : OGRSpatialReference m_oSRS{};
91 :
92 : public:
93 : NSIDCbinDataset();
94 : ~NSIDCbinDataset() override;
95 : CPLErr GetGeoTransform(double *) override;
96 :
97 : const OGRSpatialReference *GetSpatialRef() const override;
98 : static GDALDataset *Open(GDALOpenInfo *);
99 : static int Identify(GDALOpenInfo *);
100 : };
101 :
102 4 : static const char *stripLeadingSpaces_nsidc(const char *buf)
103 : {
104 4 : const char *ptr = buf;
105 : /* Go until we run out of characters or hit something non-zero */
106 9 : while (*ptr == ' ')
107 : {
108 5 : ptr++;
109 : }
110 4 : return ptr;
111 : }
112 :
113 : /************************************************************************/
114 : /* ==================================================================== */
115 : /* NSIDCbinRasterBand */
116 : /* ==================================================================== */
117 : /************************************************************************/
118 :
119 : class NSIDCbinRasterBand final : public RawRasterBand
120 : {
121 : friend class NSIDCbinDataset;
122 :
123 : CPL_DISALLOW_COPY_ASSIGN(NSIDCbinRasterBand)
124 :
125 : public:
126 : NSIDCbinRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
127 : vsi_l_offset nImgOffset, int nPixelOffset,
128 : int nLineOffset, GDALDataType eDataType);
129 : ~NSIDCbinRasterBand() override;
130 :
131 : double GetNoDataValue(int *pbSuccess = nullptr) override;
132 : double GetScale(int *pbSuccess = nullptr) override;
133 : const char *GetUnitType() override;
134 : };
135 :
136 : /************************************************************************/
137 : /* NSIDCbinRasterBand() */
138 : /************************************************************************/
139 :
140 1 : NSIDCbinRasterBand::NSIDCbinRasterBand(GDALDataset *poDSIn, int nBandIn,
141 : VSILFILE *fpRawIn,
142 : vsi_l_offset nImgOffsetIn,
143 : int nPixelOffsetIn, int nLineOffsetIn,
144 1 : GDALDataType eDataTypeIn)
145 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
146 : nLineOffsetIn, eDataTypeIn,
147 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
148 1 : RawRasterBand::OwnFP::NO)
149 : {
150 1 : }
151 :
152 : /************************************************************************/
153 : /* ~NSIDCbinRasterBand() */
154 : /************************************************************************/
155 :
156 2 : NSIDCbinRasterBand::~NSIDCbinRasterBand()
157 : {
158 2 : }
159 :
160 : /************************************************************************/
161 : /* GetNoDataValue() */
162 : /************************************************************************/
163 :
164 0 : double NSIDCbinRasterBand::GetNoDataValue(int *pbSuccess)
165 :
166 : {
167 0 : if (pbSuccess != nullptr)
168 0 : *pbSuccess = TRUE;
169 : // we might check this if other format variants can be different
170 : // or if we change the Band type, or if we generalize to choosing Byte vs.
171 : // Float type but for now it's constant
172 : // https://lists.osgeo.org/pipermail/gdal-dev/2022-August/056144.html
173 : // const char *pszLine = poPDS->sHeader.missing_int;
174 0 : return 255.0; // CPLAtof(pszLine);
175 : }
176 :
177 : /************************************************************************/
178 : /* GetScale() */
179 : /************************************************************************/
180 :
181 0 : double NSIDCbinRasterBand::GetScale(int *pbSuccess)
182 : {
183 0 : if (pbSuccess != nullptr)
184 0 : *pbSuccess = TRUE;
185 : // again just use a constant unless we see other file variants
186 : // also, this might be fraction rather than percentage
187 : // atof(reinterpret_cast<NSIDCbinDataset*>(poDS)->sHeader.scaling)/100;
188 0 : return 0.4;
189 : }
190 :
191 : /************************************************************************/
192 : /* NSIDCbinDataset() */
193 : /************************************************************************/
194 :
195 1 : NSIDCbinDataset::NSIDCbinDataset() : fp(nullptr), m_oSRS(OGRSpatialReference())
196 : {
197 1 : adfGeoTransform[0] = 0.0;
198 1 : adfGeoTransform[1] = 1.0;
199 1 : adfGeoTransform[2] = 0.0;
200 1 : adfGeoTransform[3] = 0.0;
201 1 : adfGeoTransform[4] = 0.0;
202 1 : adfGeoTransform[5] = 1.0;
203 1 : }
204 :
205 : /************************************************************************/
206 : /* ~NSIDCbinDataset() */
207 : /************************************************************************/
208 :
209 2 : NSIDCbinDataset::~NSIDCbinDataset()
210 :
211 : {
212 1 : if (fp)
213 1 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
214 1 : fp = nullptr;
215 2 : }
216 :
217 : /************************************************************************/
218 : /* Open() */
219 : /************************************************************************/
220 :
221 29910 : GDALDataset *NSIDCbinDataset::Open(GDALOpenInfo *poOpenInfo)
222 :
223 : {
224 :
225 : // Confirm that the header is compatible with a NSIDC dataset.
226 29910 : if (!Identify(poOpenInfo))
227 29880 : return nullptr;
228 :
229 : // Confirm the requested access is supported.
230 4 : if (poOpenInfo->eAccess == GA_Update)
231 : {
232 0 : ReportUpdateNotSupportedByDriver("NSIDCbin");
233 0 : return nullptr;
234 : }
235 :
236 : // Check that the file pointer from GDALOpenInfo* is available
237 4 : if (poOpenInfo->fpL == nullptr)
238 : {
239 0 : return nullptr;
240 : }
241 :
242 : /* -------------------------------------------------------------------- */
243 : /* Create a corresponding GDALDataset. */
244 : /* -------------------------------------------------------------------- */
245 5 : auto poDS = std::make_unique<NSIDCbinDataset>();
246 :
247 1 : poDS->eAccess = poOpenInfo->eAccess;
248 1 : std::swap(poDS->fp, poOpenInfo->fpL);
249 :
250 : /* -------------------------------------------------------------------- */
251 : /* Read the header information. */
252 : /* -------------------------------------------------------------------- */
253 1 : if (VSIFReadL(&(poDS->sHeader), 300, 1, poDS->fp) != 1)
254 : {
255 0 : CPLError(CE_Failure, CPLE_FileIO,
256 : "Attempt to read 300 byte header filed on file %s\n",
257 : poOpenInfo->pszFilename);
258 0 : return nullptr;
259 : }
260 :
261 : // avoid unused warnings
262 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.missing_int);
263 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.internal1);
264 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.latitude);
265 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.greenwich);
266 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.internal2);
267 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.jpole);
268 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.ipole);
269 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.julian_start);
270 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.hour_start);
271 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.minute_start);
272 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.julian_end);
273 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.hour_end);
274 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.minute_end);
275 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.channel);
276 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.scaling);
277 :
278 : /* -------------------------------------------------------------------- */
279 : /* Extract information of interest from the header. */
280 : /* -------------------------------------------------------------------- */
281 :
282 1 : poDS->nRasterXSize = atoi(poDS->sHeader.columns);
283 1 : poDS->nRasterYSize = atoi(poDS->sHeader.rows);
284 :
285 1 : const char *psHeader = reinterpret_cast<char *>(poOpenInfo->pabyHeader);
286 1 : bool south = STARTS_WITH(psHeader + 230, "ANTARCTIC");
287 :
288 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
289 : {
290 0 : return nullptr;
291 : }
292 :
293 : /* -------------------------------------------------------------------- */
294 : /* Extract metadata from the header. */
295 : /* -------------------------------------------------------------------- */
296 :
297 1 : poDS->SetMetadataItem("INSTRUMENT", poDS->sHeader.instrument);
298 1 : poDS->SetMetadataItem("YEAR", stripLeadingSpaces_nsidc(poDS->sHeader.year));
299 2 : poDS->SetMetadataItem("JULIAN_DAY",
300 1 : stripLeadingSpaces_nsidc(poDS->sHeader.julian));
301 2 : poDS->SetMetadataItem(
302 : "DATA_DESCRIPTORS",
303 1 : stripLeadingSpaces_nsidc(poDS->sHeader.data_descriptors));
304 1 : poDS->SetMetadataItem("IMAGE_TITLE", poDS->sHeader.imagetitle);
305 2 : poDS->SetMetadataItem("FILENAME",
306 1 : stripLeadingSpaces_nsidc(poDS->sHeader.filename));
307 1 : poDS->SetMetadataItem("DATA_INFORMATION", poDS->sHeader.data_information);
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Create band information objects. */
311 : /* -------------------------------------------------------------------- */
312 1 : int nBytesPerSample = 1;
313 :
314 : auto poBand = std::make_unique<NSIDCbinRasterBand>(
315 1 : poDS.get(), 1, poDS->fp, 300, nBytesPerSample, poDS->nRasterXSize,
316 3 : GDT_Byte);
317 1 : if (!poBand->IsValid())
318 0 : return nullptr;
319 1 : poDS->SetBand(1, std::move(poBand));
320 :
321 : /* -------------------------------------------------------------------- */
322 : /* Geotransform, we simply know this from the documentation. */
323 : /* If we have similar binary files (at 12.5km for example) then */
324 : /* need more nuanced handling */
325 : /* Projection, this is not technically enough, because the old */
326 : /* stuff is Hughes 1980. */
327 : /* FIXME: old or new epsg codes based on header info, or jul/year */
328 : /* -------------------------------------------------------------------- */
329 :
330 1 : int epsg = -1;
331 1 : if (south)
332 : {
333 1 : poDS->adfGeoTransform[0] = -3950000.0;
334 1 : poDS->adfGeoTransform[1] = 25000;
335 1 : poDS->adfGeoTransform[2] = 0.0;
336 1 : poDS->adfGeoTransform[3] = 4350000.0;
337 1 : poDS->adfGeoTransform[4] = 0.0;
338 1 : poDS->adfGeoTransform[5] = -25000;
339 :
340 1 : epsg = 3976;
341 : }
342 : else
343 : {
344 0 : poDS->adfGeoTransform[0] = -3837500;
345 0 : poDS->adfGeoTransform[1] = 25000;
346 0 : poDS->adfGeoTransform[2] = 0.0;
347 0 : poDS->adfGeoTransform[3] = 5837500;
348 0 : poDS->adfGeoTransform[4] = 0.0;
349 0 : poDS->adfGeoTransform[5] = -25000;
350 :
351 0 : epsg = 3413;
352 : }
353 :
354 1 : if (poDS->m_oSRS.importFromEPSG(epsg) != OGRERR_NONE)
355 : {
356 0 : CPLError(CE_Failure, CPLE_AppDefined,
357 : "Unknown error initializing SRS from ESPG code. ");
358 0 : return nullptr;
359 : }
360 :
361 : /* -------------------------------------------------------------------- */
362 : /* Initialize any PAM information. */
363 : /* -------------------------------------------------------------------- */
364 1 : poDS->SetDescription(poOpenInfo->pszFilename);
365 1 : poDS->TryLoadXML();
366 :
367 1 : return poDS.release();
368 : }
369 :
370 : /************************************************************************/
371 : /* Identify() */
372 : /************************************************************************/
373 29914 : int NSIDCbinDataset::Identify(GDALOpenInfo *poOpenInfo)
374 : {
375 :
376 : // -------------------------------------------------------------------- /
377 : // Works for daily and monthly, north and south NSIDC binary files /
378 : // north and south are different dimensions, different extents but /
379 : // both are 25000m resolution.
380 : //
381 : // First we check to see if the file has the expected header /
382 : // bytes. /
383 : // -------------------------------------------------------------------- /
384 :
385 29914 : if (poOpenInfo->nHeaderBytes < 300 || poOpenInfo->fpL == nullptr)
386 28034 : return FALSE;
387 :
388 1880 : const char *psHeader = reinterpret_cast<char *>(poOpenInfo->pabyHeader);
389 : // Check if century values seem reasonable.
390 1880 : if (!(EQUALN(psHeader + 103, "20", 2) || EQUALN(psHeader + 103, "19", 2) ||
391 : // the first files 1978 don't have a space at the start
392 1879 : EQUALN(psHeader + 102, "20", 2) || EQUALN(psHeader + 102, "19", 2)))
393 : {
394 1879 : return FALSE;
395 : }
396 :
397 : // Check if descriptors reasonable.
398 1 : if (!(STARTS_WITH(psHeader + 230, "ANTARCTIC") ||
399 0 : STARTS_WITH(psHeader + 230, "ARCTIC")))
400 : {
401 :
402 0 : return FALSE;
403 : }
404 :
405 1 : return TRUE;
406 : }
407 :
408 : /************************************************************************/
409 : /* GetSpatialRef() */
410 : /************************************************************************/
411 :
412 0 : const OGRSpatialReference *NSIDCbinDataset::GetSpatialRef() const
413 : {
414 0 : return &m_oSRS;
415 : }
416 :
417 : /************************************************************************/
418 : /* GetGeoTransform() */
419 : /************************************************************************/
420 :
421 0 : CPLErr NSIDCbinDataset::GetGeoTransform(double *padfTransform)
422 :
423 : {
424 0 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
425 :
426 0 : return CE_None;
427 : }
428 :
429 : /************************************************************************/
430 : /* GetUnitType() */
431 : /************************************************************************/
432 :
433 0 : const char *NSIDCbinRasterBand::GetUnitType()
434 : {
435 : // undecided, atm stick with Byte but may switch to Float and lose values >
436 : // 250 or generalize to non-raw driver
437 : // https://lists.osgeo.org/pipermail/gdal-dev/2022-August/056144.html
438 : // if (eDataType == GDT_Float32)
439 : // return "Percentage"; // or "Fraction [0,1]"
440 :
441 : // Byte values don't have a clear unit type
442 0 : return "";
443 : }
444 :
445 : /************************************************************************/
446 : /* GDALRegister_NSIDCbin() */
447 : /************************************************************************/
448 :
449 1686 : void GDALRegister_NSIDCbin()
450 :
451 : {
452 1686 : if (GDALGetDriverByName("NSIDCbin") != nullptr)
453 302 : return;
454 :
455 1384 : GDALDriver *poDriver = new GDALDriver();
456 :
457 1384 : poDriver->SetDescription("NSIDCbin");
458 1384 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
459 1384 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
460 1384 : "NSIDC Sea Ice Concentrations binary (.bin)");
461 1384 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
462 1384 : "drivers/raster/nsidcbin.html");
463 1384 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
464 1384 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bin");
465 :
466 1384 : poDriver->pfnOpen = NSIDCbinDataset::Open;
467 :
468 1384 : GetGDALDriverManager()->RegisterDriver(poDriver);
469 : }
|
