Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  NDF Driver
       4             :  * Purpose:  Implementation of NLAPS Data Format read support.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2005, Frank Warmerdam
       9             :  * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_string.h"
      15             : #include "gdal_frmts.h"
      16             : #include "ogr_spatialref.h"
      17             : #include "rawdataset.h"
      18             : 
      19             : /************************************************************************/
      20             : /* ==================================================================== */
      21             : /*                              NDFDataset                              */
      22             : /* ==================================================================== */
      23             : /************************************************************************/
      24             : 
      25             : class NDFDataset final : public RawDataset
      26             : {
      27             :     double adfGeoTransform[6];
      28             : 
      29             :     OGRSpatialReference m_oSRS{};
      30             :     char **papszExtraFiles;
      31             : 
      32             :     char **papszHeader;
      33             :     const char *Get(const char *pszKey, const char *pszDefault);
      34             : 
      35             :     CPL_DISALLOW_COPY_ASSIGN(NDFDataset)
      36             : 
      37             :     CPLErr Close() override;
      38             : 
      39             :   public:
      40             :     NDFDataset();
      41             :     ~NDFDataset() override;
      42             : 
      43             :     CPLErr GetGeoTransform(double *padfTransform) override;
      44             : 
      45           1 :     const OGRSpatialReference *GetSpatialRef() const override
      46             :     {
      47           1 :         return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
      48             :     }
      49             : 
      50             :     char **GetFileList(void) override;
      51             : 
      52             :     static GDALDataset *Open(GDALOpenInfo *);
      53             :     static int Identify(GDALOpenInfo *);
      54             : };
      55             : 
      56             : /************************************************************************/
      57             : /*                            NDFDataset()                             */
      58             : /************************************************************************/
      59             : 
      60           2 : NDFDataset::NDFDataset() : papszExtraFiles(nullptr), papszHeader(nullptr)
      61             : {
      62           2 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
      63           2 :     adfGeoTransform[0] = 0.0;
      64           2 :     adfGeoTransform[1] = 1.0;
      65           2 :     adfGeoTransform[2] = 0.0;
      66           2 :     adfGeoTransform[3] = 0.0;
      67           2 :     adfGeoTransform[4] = 0.0;
      68           2 :     adfGeoTransform[5] = 1.0;
      69           2 : }
      70             : 
      71             : /************************************************************************/
      72             : /*                            ~NDFDataset()                             */
      73             : /************************************************************************/
      74             : 
      75           4 : NDFDataset::~NDFDataset()
      76             : 
      77             : {
      78           2 :     NDFDataset::Close();
      79           4 : }
      80             : 
      81             : /************************************************************************/
      82             : /*                              Close()                                 */
      83             : /************************************************************************/
      84             : 
      85           4 : CPLErr NDFDataset::Close()
      86             : {
      87           4 :     CPLErr eErr = CE_None;
      88           4 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
      89             :     {
      90           2 :         if (NDFDataset::FlushCache(true) != CE_None)
      91           0 :             eErr = CE_Failure;
      92             : 
      93           2 :         CSLDestroy(papszHeader);
      94           2 :         CSLDestroy(papszExtraFiles);
      95             : 
      96           2 :         if (GDALPamDataset::Close() != CE_None)
      97           0 :             eErr = CE_Failure;
      98             :     }
      99           4 :     return eErr;
     100             : }
     101             : 
     102             : /************************************************************************/
     103             : /*                          GetGeoTransform()                           */
     104             : /************************************************************************/
     105             : 
     106           1 : CPLErr NDFDataset::GetGeoTransform(double *padfTransform)
     107             : 
     108             : {
     109           1 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     110           1 :     return CE_None;
     111             : }
     112             : 
     113             : /************************************************************************/
     114             : /*                                Get()                                 */
     115             : /*                                                                      */
     116             : /*      Fetch a value from the header by keyword.                       */
     117             : /************************************************************************/
     118             : 
     119          26 : const char *NDFDataset::Get(const char *pszKey, const char *pszDefault)
     120             : 
     121             : {
     122          26 :     const char *pszResult = CSLFetchNameValue(papszHeader, pszKey);
     123             : 
     124          26 :     if (pszResult == nullptr)
     125           0 :         return pszDefault;
     126             : 
     127          26 :     return pszResult;
     128             : }
     129             : 
     130             : /************************************************************************/
     131             : /*                            GetFileList()                             */
     132             : /************************************************************************/
     133             : 
     134           1 : char **NDFDataset::GetFileList()
     135             : 
     136             : {
     137             :     // Main data file, etc.
     138           1 :     char **papszFileList = GDALPamDataset::GetFileList();
     139             : 
     140             :     // Header file.
     141           1 :     papszFileList = CSLInsertStrings(papszFileList, -1, papszExtraFiles);
     142             : 
     143           1 :     return papszFileList;
     144             : }
     145             : 
     146             : /************************************************************************/
     147             : /*                            Identify()                                */
     148             : /************************************************************************/
     149             : 
     150       51572 : int NDFDataset::Identify(GDALOpenInfo *poOpenInfo)
     151             : 
     152             : {
     153             :     /* -------------------------------------------------------------------- */
     154             :     /*      The user must select the header file (i.e. .H1).                */
     155             :     /* -------------------------------------------------------------------- */
     156       55073 :     return poOpenInfo->nHeaderBytes >= 50 &&
     157        3501 :            (STARTS_WITH_CI(
     158             :                 reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
     159        3497 :                 "NDF_REVISION=2") ||
     160        3497 :             STARTS_WITH_CI(
     161             :                 reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
     162       51572 :                 "NDF_REVISION=0"));
     163             : }
     164             : 
     165             : /************************************************************************/
     166             : /*                                Open()                                */
     167             : /************************************************************************/
     168             : 
     169           2 : GDALDataset *NDFDataset::Open(GDALOpenInfo *poOpenInfo)
     170             : 
     171             : {
     172             :     /* -------------------------------------------------------------------- */
     173             :     /*      The user must select the header file (i.e. .H1).                */
     174             :     /* -------------------------------------------------------------------- */
     175           2 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     176           0 :         return nullptr;
     177             : 
     178             :     /* -------------------------------------------------------------------- */
     179             :     /*      Confirm the requested access is supported.                      */
     180             :     /* -------------------------------------------------------------------- */
     181           2 :     if (poOpenInfo->eAccess == GA_Update)
     182             :     {
     183           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     184             :                  "The NDF driver does not support update access to existing"
     185             :                  " datasets.");
     186           0 :         return nullptr;
     187             :     }
     188             :     /* -------------------------------------------------------------------- */
     189             :     /*      Read and process the header into a local name/value             */
     190             :     /*      stringlist.  We just take off the trailing semicolon.  The      */
     191             :     /*      keyword is already separated from the value by an equal         */
     192             :     /*      sign.                                                           */
     193             :     /* -------------------------------------------------------------------- */
     194             : 
     195           2 :     const char *pszLine = nullptr;
     196           2 :     const int nHeaderMax = 1000;
     197           2 :     int nHeaderLines = 0;
     198             :     char **papszHeader =
     199           2 :         static_cast<char **>(CPLMalloc(sizeof(char *) * (nHeaderMax + 1)));
     200             : 
     201         210 :     while (nHeaderLines < nHeaderMax &&
     202         212 :            (pszLine = CPLReadLineL(poOpenInfo->fpL)) != nullptr &&
     203         106 :            !EQUAL(pszLine, "END_OF_HDR;"))
     204             :     {
     205         104 :         if (strstr(pszLine, "=") == nullptr)
     206           0 :             break;
     207             : 
     208         104 :         char *pszFixed = CPLStrdup(pszLine);
     209         104 :         if (pszFixed[strlen(pszFixed) - 1] == ';')
     210         104 :             pszFixed[strlen(pszFixed) - 1] = '\0';
     211             : 
     212         104 :         papszHeader[nHeaderLines++] = pszFixed;
     213         104 :         papszHeader[nHeaderLines] = nullptr;
     214             :     }
     215           2 :     CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
     216           2 :     poOpenInfo->fpL = nullptr;
     217             : 
     218           2 :     if (CSLFetchNameValue(papszHeader, "PIXELS_PER_LINE") == nullptr ||
     219           2 :         CSLFetchNameValue(papszHeader, "LINES_PER_DATA_FILE") == nullptr ||
     220           6 :         CSLFetchNameValue(papszHeader, "BITS_PER_PIXEL") == nullptr ||
     221           2 :         CSLFetchNameValue(papszHeader, "PIXEL_FORMAT") == nullptr)
     222             :     {
     223           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     224             :                  "Dataset appears to be NDF but is missing a required field.");
     225           0 :         CSLDestroy(papszHeader);
     226           0 :         return nullptr;
     227             :     }
     228             : 
     229           4 :     if (!EQUAL(CSLFetchNameValue(papszHeader, "PIXEL_FORMAT"), "BYTE") ||
     230           2 :         !EQUAL(CSLFetchNameValue(papszHeader, "BITS_PER_PIXEL"), "8"))
     231             :     {
     232           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     233             :                  "Currently NDF driver supports only 8bit BYTE format.");
     234           0 :         CSLDestroy(papszHeader);
     235           0 :         return nullptr;
     236             :     }
     237             : 
     238             :     /* -------------------------------------------------------------------- */
     239             :     /*      Confirm the requested access is supported.                      */
     240             :     /* -------------------------------------------------------------------- */
     241           2 :     if (poOpenInfo->eAccess == GA_Update)
     242             :     {
     243           0 :         CSLDestroy(papszHeader);
     244           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     245             :                  "The NDF driver does not support update access to existing"
     246             :                  " datasets.\n");
     247           0 :         return nullptr;
     248             :     }
     249             : 
     250             :     /* -------------------------------------------------------------------- */
     251             :     /*      Create a corresponding GDALDataset.                             */
     252             :     /* -------------------------------------------------------------------- */
     253           4 :     auto poDS = std::make_unique<NDFDataset>();
     254           2 :     poDS->papszHeader = papszHeader;
     255             : 
     256           2 :     poDS->nRasterXSize = atoi(poDS->Get("PIXELS_PER_LINE", ""));
     257           2 :     poDS->nRasterYSize = atoi(poDS->Get("LINES_PER_DATA_FILE", ""));
     258             : 
     259             :     /* -------------------------------------------------------------------- */
     260             :     /*      Create a raw raster band for each file.                         */
     261             :     /* -------------------------------------------------------------------- */
     262             :     const char *pszBand =
     263           2 :         CSLFetchNameValue(papszHeader, "NUMBER_OF_BANDS_IN_VOLUME");
     264           2 :     if (pszBand == nullptr)
     265             :     {
     266           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot find band count");
     267           0 :         return nullptr;
     268             :     }
     269           2 :     const int nBands = atoi(pszBand);
     270             : 
     271           4 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
     272           2 :         !GDALCheckBandCount(nBands, FALSE))
     273             :     {
     274           0 :         return nullptr;
     275             :     }
     276             : 
     277           4 :     for (int iBand = 0; iBand < nBands; iBand++)
     278             :     {
     279             :         char szKey[100];
     280           2 :         snprintf(szKey, sizeof(szKey), "BAND%d_FILENAME", iBand + 1);
     281           2 :         CPLString osFilename = poDS->Get(szKey, "");
     282             : 
     283             :         // NDF1 file do not include the band filenames.
     284           2 :         if (osFilename.empty())
     285             :         {
     286             :             char szBandExtension[15];
     287           0 :             snprintf(szBandExtension, sizeof(szBandExtension), "I%d",
     288             :                      iBand + 1);
     289             :             osFilename =
     290           0 :                 CPLResetExtension(poOpenInfo->pszFilename, szBandExtension);
     291             :         }
     292             :         else
     293             :         {
     294           2 :             CPLString osBasePath = CPLGetPath(poOpenInfo->pszFilename);
     295           2 :             osFilename = CPLFormFilename(osBasePath, osFilename, nullptr);
     296             :         }
     297             : 
     298           2 :         VSILFILE *fpRaw = VSIFOpenL(osFilename, "rb");
     299           2 :         if (fpRaw == nullptr)
     300             :         {
     301           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     302             :                      "Failed to open band file: %s", osFilename.c_str());
     303           0 :             return nullptr;
     304             :         }
     305           2 :         poDS->papszExtraFiles = CSLAddString(poDS->papszExtraFiles, osFilename);
     306             : 
     307             :         auto poBand = RawRasterBand::Create(
     308           2 :             poDS.get(), iBand + 1, fpRaw, 0, 1, poDS->nRasterXSize, GDT_Byte,
     309             :             RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     310           2 :             RawRasterBand::OwnFP::YES);
     311           2 :         if (!poBand)
     312           0 :             return nullptr;
     313             : 
     314           2 :         snprintf(szKey, sizeof(szKey), "BAND%d_NAME", iBand + 1);
     315           2 :         poBand->SetDescription(poDS->Get(szKey, ""));
     316             : 
     317           2 :         snprintf(szKey, sizeof(szKey), "BAND%d_WAVELENGTHS", iBand + 1);
     318           2 :         poBand->SetMetadataItem("WAVELENGTHS", poDS->Get(szKey, ""));
     319             : 
     320           2 :         snprintf(szKey, sizeof(szKey), "BAND%d_RADIOMETRIC_GAINS/BIAS",
     321             :                  iBand + 1);
     322           2 :         poBand->SetMetadataItem("RADIOMETRIC_GAINS_BIAS", poDS->Get(szKey, ""));
     323             : 
     324           2 :         poDS->SetBand(iBand + 1, std::move(poBand));
     325             :     }
     326             : 
     327             :     /* -------------------------------------------------------------------- */
     328             :     /*      Fetch and parse USGS projection parameters.                     */
     329             :     /* -------------------------------------------------------------------- */
     330           2 :     double adfUSGSParams[15] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
     331             :     const CPLStringList aosParamTokens(CSLTokenizeStringComplex(
     332           4 :         poDS->Get("USGS_PROJECTION_NUMBER", ""), ",", FALSE, TRUE));
     333             : 
     334           2 :     if (aosParamTokens.size() >= 15)
     335             :     {
     336           0 :         for (int i = 0; i < 15; i++)
     337           0 :             adfUSGSParams[i] = CPLAtof(aosParamTokens[i]);
     338             :     }
     339             : 
     340             :     /* -------------------------------------------------------------------- */
     341             :     /*      Minimal georef support ... should add full USGS style           */
     342             :     /*      support at some point.                                          */
     343             :     /* -------------------------------------------------------------------- */
     344           2 :     const int nUSGSProjection = atoi(poDS->Get("USGS_PROJECTION_NUMBER", ""));
     345           2 :     const int nZone = atoi(poDS->Get("USGS_MAP_ZONE", "0"));
     346             : 
     347           4 :     OGRSpatialReference oSRS;
     348           2 :     oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     349           2 :     oSRS.importFromUSGS(nUSGSProjection, nZone, adfUSGSParams, 12);
     350             : 
     351           4 :     const CPLString osDatum = poDS->Get("HORIZONTAL_DATUM", "");
     352           2 :     if (EQUAL(osDatum, "WGS84") || EQUAL(osDatum, "NAD83") ||
     353           0 :         EQUAL(osDatum, "NAD27"))
     354             :     {
     355           2 :         oSRS.SetWellKnownGeogCS(osDatum);
     356             :     }
     357           0 :     else if (STARTS_WITH_CI(osDatum, "NAD27"))
     358             :     {
     359           0 :         oSRS.SetWellKnownGeogCS("NAD27");
     360             :     }
     361             :     else
     362             :     {
     363           0 :         CPLError(CE_Warning, CPLE_AppDefined,
     364             :                  "Unrecognized datum name in NLAPS/NDF file:%s, "
     365             :                  "assuming WGS84.",
     366             :                  osDatum.c_str());
     367           0 :         oSRS.SetWellKnownGeogCS("WGS84");
     368             :     }
     369             : 
     370           2 :     if (oSRS.GetRoot() != nullptr)
     371             :     {
     372           2 :         poDS->m_oSRS = std::move(oSRS);
     373             :     }
     374             : 
     375             :     /* -------------------------------------------------------------------- */
     376             :     /*      Get geotransform.                                               */
     377             :     /* -------------------------------------------------------------------- */
     378             :     char **papszUL =
     379           2 :         CSLTokenizeString2(poDS->Get("UPPER_LEFT_CORNER", ""), ",", 0);
     380             :     char **papszUR =
     381           2 :         CSLTokenizeString2(poDS->Get("UPPER_RIGHT_CORNER", ""), ",", 0);
     382             :     char **papszLL =
     383           2 :         CSLTokenizeString2(poDS->Get("LOWER_LEFT_CORNER", ""), ",", 0);
     384             : 
     385           4 :     if (CSLCount(papszUL) == 4 && CSLCount(papszUR) == 4 &&
     386           2 :         CSLCount(papszLL) == 4)
     387             :     {
     388           2 :         poDS->adfGeoTransform[0] = CPLAtof(papszUL[2]);
     389           2 :         poDS->adfGeoTransform[1] = (CPLAtof(papszUR[2]) - CPLAtof(papszUL[2])) /
     390           2 :                                    (poDS->nRasterXSize - 1);
     391           2 :         poDS->adfGeoTransform[2] = (CPLAtof(papszUR[3]) - CPLAtof(papszUL[3])) /
     392           2 :                                    (poDS->nRasterXSize - 1);
     393             : 
     394           2 :         poDS->adfGeoTransform[3] = CPLAtof(papszUL[3]);
     395           2 :         poDS->adfGeoTransform[4] = (CPLAtof(papszLL[2]) - CPLAtof(papszUL[2])) /
     396           2 :                                    (poDS->nRasterYSize - 1);
     397           2 :         poDS->adfGeoTransform[5] = (CPLAtof(papszLL[3]) - CPLAtof(papszUL[3])) /
     398           2 :                                    (poDS->nRasterYSize - 1);
     399             : 
     400             :         // Move origin up-left half a pixel.
     401           2 :         poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
     402           2 :         poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[4] * 0.5;
     403           2 :         poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[2] * 0.5;
     404           2 :         poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[5] * 0.5;
     405             :     }
     406             : 
     407           2 :     CSLDestroy(papszUL);
     408           2 :     CSLDestroy(papszLL);
     409           2 :     CSLDestroy(papszUR);
     410             : 
     411             :     /* -------------------------------------------------------------------- */
     412             :     /*      Initialize any PAM information.                                 */
     413             :     /* -------------------------------------------------------------------- */
     414           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     415           2 :     poDS->TryLoadXML();
     416             : 
     417             :     /* -------------------------------------------------------------------- */
     418             :     /*      Check for overviews.                                            */
     419             :     /* -------------------------------------------------------------------- */
     420           2 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     421             : 
     422           2 :     return poDS.release();
     423             : }
     424             : 
     425             : /************************************************************************/
     426             : /*                          GDALRegister_NDF()                          */
     427             : /************************************************************************/
     428             : 
     429        1595 : void GDALRegister_NDF()
     430             : 
     431             : {
     432        1595 :     if (GDALGetDriverByName("NDF") != nullptr)
     433         302 :         return;
     434             : 
     435        1293 :     GDALDriver *poDriver = new GDALDriver();
     436             : 
     437        1293 :     poDriver->SetDescription("NDF");
     438        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     439        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NLAPS Data Format");
     440        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ndf.html");
     441        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     442             : 
     443        1293 :     poDriver->pfnIdentify = NDFDataset::Identify;
     444        1293 :     poDriver->pfnOpen = NDFDataset::Open;
     445             : 
     446        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     447             : }