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