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