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