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 : GDALGeoTransform m_gt{};
89 : CPL_DISALLOW_COPY_ASSIGN(NSIDCbinDataset)
90 : OGRSpatialReference m_oSRS{};
91 :
92 : public:
93 : NSIDCbinDataset();
94 : ~NSIDCbinDataset() override;
95 : CPLErr GetGeoTransform(GDALGeoTransform >) const 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 : }
198 :
199 : /************************************************************************/
200 : /* ~NSIDCbinDataset() */
201 : /************************************************************************/
202 :
203 2 : NSIDCbinDataset::~NSIDCbinDataset()
204 :
205 : {
206 1 : if (fp)
207 1 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
208 1 : fp = nullptr;
209 2 : }
210 :
211 : /************************************************************************/
212 : /* Open() */
213 : /************************************************************************/
214 :
215 32306 : GDALDataset *NSIDCbinDataset::Open(GDALOpenInfo *poOpenInfo)
216 :
217 : {
218 :
219 : // Confirm that the header is compatible with a NSIDC dataset.
220 32306 : if (!Identify(poOpenInfo))
221 32295 : return nullptr;
222 :
223 : // Confirm the requested access is supported.
224 9 : if (poOpenInfo->eAccess == GA_Update)
225 : {
226 0 : ReportUpdateNotSupportedByDriver("NSIDCbin");
227 0 : return nullptr;
228 : }
229 :
230 : // Check that the file pointer from GDALOpenInfo* is available
231 9 : if (poOpenInfo->fpL == nullptr)
232 : {
233 0 : return nullptr;
234 : }
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Create a corresponding GDALDataset. */
238 : /* -------------------------------------------------------------------- */
239 10 : auto poDS = std::make_unique<NSIDCbinDataset>();
240 :
241 1 : poDS->eAccess = poOpenInfo->eAccess;
242 1 : std::swap(poDS->fp, poOpenInfo->fpL);
243 :
244 : /* -------------------------------------------------------------------- */
245 : /* Read the header information. */
246 : /* -------------------------------------------------------------------- */
247 1 : if (VSIFReadL(&(poDS->sHeader), 300, 1, poDS->fp) != 1)
248 : {
249 0 : CPLError(CE_Failure, CPLE_FileIO,
250 : "Attempt to read 300 byte header filed on file %s\n",
251 : poOpenInfo->pszFilename);
252 0 : return nullptr;
253 : }
254 :
255 : // avoid unused warnings
256 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.missing_int);
257 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.internal1);
258 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.latitude);
259 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.greenwich);
260 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.internal2);
261 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.jpole);
262 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.ipole);
263 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.julian_start);
264 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.hour_start);
265 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.minute_start);
266 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.julian_end);
267 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.hour_end);
268 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.minute_end);
269 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.channel);
270 1 : CPL_IGNORE_RET_VAL(poDS->sHeader.scaling);
271 :
272 : /* -------------------------------------------------------------------- */
273 : /* Extract information of interest from the header. */
274 : /* -------------------------------------------------------------------- */
275 :
276 1 : poDS->nRasterXSize = atoi(poDS->sHeader.columns);
277 1 : poDS->nRasterYSize = atoi(poDS->sHeader.rows);
278 :
279 1 : const char *psHeader = reinterpret_cast<char *>(poOpenInfo->pabyHeader);
280 1 : bool south = STARTS_WITH(psHeader + 230, "ANTARCTIC");
281 :
282 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
283 : {
284 0 : return nullptr;
285 : }
286 :
287 : /* -------------------------------------------------------------------- */
288 : /* Extract metadata from the header. */
289 : /* -------------------------------------------------------------------- */
290 :
291 1 : poDS->SetMetadataItem("INSTRUMENT", poDS->sHeader.instrument);
292 1 : poDS->SetMetadataItem("YEAR", stripLeadingSpaces_nsidc(poDS->sHeader.year));
293 2 : poDS->SetMetadataItem("JULIAN_DAY",
294 1 : stripLeadingSpaces_nsidc(poDS->sHeader.julian));
295 2 : poDS->SetMetadataItem(
296 : "DATA_DESCRIPTORS",
297 1 : stripLeadingSpaces_nsidc(poDS->sHeader.data_descriptors));
298 1 : poDS->SetMetadataItem("IMAGE_TITLE", poDS->sHeader.imagetitle);
299 2 : poDS->SetMetadataItem("FILENAME",
300 1 : stripLeadingSpaces_nsidc(poDS->sHeader.filename));
301 1 : poDS->SetMetadataItem("DATA_INFORMATION", poDS->sHeader.data_information);
302 :
303 : /* -------------------------------------------------------------------- */
304 : /* Create band information objects. */
305 : /* -------------------------------------------------------------------- */
306 1 : int nBytesPerSample = 1;
307 :
308 : auto poBand = std::make_unique<NSIDCbinRasterBand>(
309 1 : poDS.get(), 1, poDS->fp, 300, nBytesPerSample, poDS->nRasterXSize,
310 3 : GDT_Byte);
311 1 : if (!poBand->IsValid())
312 0 : return nullptr;
313 1 : poDS->SetBand(1, std::move(poBand));
314 :
315 : /* -------------------------------------------------------------------- */
316 : /* Geotransform, we simply know this from the documentation. */
317 : /* If we have similar binary files (at 12.5km for example) then */
318 : /* need more nuanced handling */
319 : /* Projection, this is not technically enough, because the old */
320 : /* stuff is Hughes 1980. */
321 : /* FIXME: old or new epsg codes based on header info, or jul/year */
322 : /* -------------------------------------------------------------------- */
323 :
324 1 : int epsg = -1;
325 1 : if (south)
326 : {
327 1 : poDS->m_gt[0] = -3950000.0;
328 1 : poDS->m_gt[1] = 25000;
329 1 : poDS->m_gt[2] = 0.0;
330 1 : poDS->m_gt[3] = 4350000.0;
331 1 : poDS->m_gt[4] = 0.0;
332 1 : poDS->m_gt[5] = -25000;
333 :
334 1 : epsg = 3976;
335 : }
336 : else
337 : {
338 0 : poDS->m_gt[0] = -3837500;
339 0 : poDS->m_gt[1] = 25000;
340 0 : poDS->m_gt[2] = 0.0;
341 0 : poDS->m_gt[3] = 5837500;
342 0 : poDS->m_gt[4] = 0.0;
343 0 : poDS->m_gt[5] = -25000;
344 :
345 0 : epsg = 3413;
346 : }
347 :
348 1 : if (poDS->m_oSRS.importFromEPSG(epsg) != OGRERR_NONE)
349 : {
350 0 : CPLError(CE_Failure, CPLE_AppDefined,
351 : "Unknown error initializing SRS from ESPG code. ");
352 0 : return nullptr;
353 : }
354 :
355 : /* -------------------------------------------------------------------- */
356 : /* Initialize any PAM information. */
357 : /* -------------------------------------------------------------------- */
358 1 : poDS->SetDescription(poOpenInfo->pszFilename);
359 1 : poDS->TryLoadXML();
360 :
361 1 : return poDS.release();
362 : }
363 :
364 : /************************************************************************/
365 : /* Identify() */
366 : /************************************************************************/
367 32306 : int NSIDCbinDataset::Identify(GDALOpenInfo *poOpenInfo)
368 : {
369 :
370 : // -------------------------------------------------------------------- /
371 : // Works for daily and monthly, north and south NSIDC binary files /
372 : // north and south are different dimensions, different extents but /
373 : // both are 25000m resolution.
374 : //
375 : // First we check to see if the file has the expected header /
376 : // bytes. /
377 : // -------------------------------------------------------------------- /
378 :
379 32306 : if (poOpenInfo->nHeaderBytes < 300 || poOpenInfo->fpL == nullptr)
380 30347 : return FALSE;
381 :
382 1959 : const char *psHeader = reinterpret_cast<char *>(poOpenInfo->pabyHeader);
383 : // Check if century values seem reasonable.
384 1959 : if (!(EQUALN(psHeader + 103, "20", 2) || EQUALN(psHeader + 103, "19", 2) ||
385 : // the first files 1978 don't have a space at the start
386 1958 : EQUALN(psHeader + 102, "20", 2) || EQUALN(psHeader + 102, "19", 2)))
387 : {
388 1958 : return FALSE;
389 : }
390 :
391 : // Check if descriptors reasonable.
392 1 : if (!(STARTS_WITH(psHeader + 230, "ANTARCTIC") ||
393 0 : STARTS_WITH(psHeader + 230, "ARCTIC")))
394 : {
395 :
396 0 : return FALSE;
397 : }
398 :
399 1 : return TRUE;
400 : }
401 :
402 : /************************************************************************/
403 : /* GetSpatialRef() */
404 : /************************************************************************/
405 :
406 0 : const OGRSpatialReference *NSIDCbinDataset::GetSpatialRef() const
407 : {
408 0 : return &m_oSRS;
409 : }
410 :
411 : /************************************************************************/
412 : /* GetGeoTransform() */
413 : /************************************************************************/
414 :
415 0 : CPLErr NSIDCbinDataset::GetGeoTransform(GDALGeoTransform >) const
416 :
417 : {
418 0 : gt = m_gt;
419 0 : return CE_None;
420 : }
421 :
422 : /************************************************************************/
423 : /* GetUnitType() */
424 : /************************************************************************/
425 :
426 0 : const char *NSIDCbinRasterBand::GetUnitType()
427 : {
428 : // undecided, atm stick with Byte but may switch to Float and lose values >
429 : // 250 or generalize to non-raw driver
430 : // https://lists.osgeo.org/pipermail/gdal-dev/2022-August/056144.html
431 : // if (eDataType == GDT_Float32)
432 : // return "Percentage"; // or "Fraction [0,1]"
433 :
434 : // Byte values don't have a clear unit type
435 0 : return "";
436 : }
437 :
438 : /************************************************************************/
439 : /* GDALRegister_NSIDCbin() */
440 : /************************************************************************/
441 :
442 1935 : void GDALRegister_NSIDCbin()
443 :
444 : {
445 1935 : if (GDALGetDriverByName("NSIDCbin") != nullptr)
446 282 : return;
447 :
448 1653 : GDALDriver *poDriver = new GDALDriver();
449 :
450 1653 : poDriver->SetDescription("NSIDCbin");
451 1653 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
452 1653 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
453 1653 : "NSIDC Sea Ice Concentrations binary (.bin)");
454 1653 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
455 1653 : "drivers/raster/nsidcbin.html");
456 1653 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
457 1653 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bin");
458 :
459 1653 : poDriver->pfnOpen = NSIDCbinDataset::Open;
460 :
461 1653 : GetGDALDriverManager()->RegisterDriver(poDriver);
462 : }
|
