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