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