Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  USGS DOQ Driver (Second Generation Format)
       4             :  * Purpose:  Implementation of DOQ2Dataset
       5             :  * Author:   Derrick J Brashear, shadow@dementia.org
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2000, Derrick J Brashear
       9             :  * Copyright (c) 2009-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 "rawdataset.h"
      33             : 
      34             : #ifndef UTM_FORMAT_defined
      35             : #define UTM_FORMAT_defined
      36             : 
      37             : static const char UTM_FORMAT[] =
      38             :     "PROJCS[\"%s / UTM zone %dN\",GEOGCS[%s,PRIMEM[\"Greenwich\",0],"
      39             :     "UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],"
      40             :     "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%d],"
      41             :     "PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],"
      42             :     "PARAMETER[\"false_northing\",0],%s]";
      43             : 
      44             : static const char WGS84_DATUM[] =
      45             :     "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]";
      46             : 
      47             : static const char WGS72_DATUM[] =
      48             :     "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]";
      49             : 
      50             : static const char NAD27_DATUM[] =
      51             :     "\"NAD27\",DATUM[\"North_American_Datum_1927\","
      52             :     "SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]";
      53             : 
      54             : static const char NAD83_DATUM[] =
      55             :     "\"NAD83\",DATUM[\"North_American_Datum_1983\","
      56             :     "SPHEROID[\"GRS 1980\",6378137,298.257222101]]";
      57             : 
      58             : #endif
      59             : 
      60             : /************************************************************************/
      61             : /* ==================================================================== */
      62             : /*                              DOQ2Dataset                             */
      63             : /* ==================================================================== */
      64             : /************************************************************************/
      65             : 
      66             : class DOQ2Dataset final : public RawDataset
      67             : {
      68             :     VSILFILE *fpImage = nullptr;  // Image data file.
      69             : 
      70             :     double dfULX = 0;
      71             :     double dfULY = 0;
      72             :     double dfXPixelSize = 0;
      73             :     double dfYPixelSize = 0;
      74             : 
      75             :     OGRSpatialReference m_oSRS{};
      76             : 
      77             :     CPL_DISALLOW_COPY_ASSIGN(DOQ2Dataset)
      78             : 
      79             :     CPLErr Close() override;
      80             : 
      81             :   public:
      82             :     DOQ2Dataset();
      83             :     ~DOQ2Dataset();
      84             : 
      85             :     CPLErr GetGeoTransform(double *padfTransform) override;
      86             : 
      87           0 :     const OGRSpatialReference *GetSpatialRef() const override
      88             :     {
      89           0 :         return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
      90             :     }
      91             : 
      92             :     static GDALDataset *Open(GDALOpenInfo *);
      93             : };
      94             : 
      95             : /************************************************************************/
      96             : /*                            DOQ2Dataset()                             */
      97             : /************************************************************************/
      98             : 
      99           1 : DOQ2Dataset::DOQ2Dataset()
     100             : {
     101           1 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     102           1 : }
     103             : 
     104             : /************************************************************************/
     105             : /*                            ~DOQ2Dataset()                            */
     106             : /************************************************************************/
     107             : 
     108           2 : DOQ2Dataset::~DOQ2Dataset()
     109             : 
     110             : {
     111           1 :     DOQ2Dataset::Close();
     112           2 : }
     113             : 
     114             : /************************************************************************/
     115             : /*                              Close()                                 */
     116             : /************************************************************************/
     117             : 
     118           2 : CPLErr DOQ2Dataset::Close()
     119             : {
     120           2 :     CPLErr eErr = CE_None;
     121           2 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     122             :     {
     123           1 :         if (DOQ2Dataset::FlushCache(true) != CE_None)
     124           0 :             eErr = CE_Failure;
     125             : 
     126           1 :         if (fpImage)
     127             :         {
     128           1 :             if (VSIFCloseL(fpImage) != 0)
     129             :             {
     130           0 :                 CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     131           0 :                 eErr = CE_Failure;
     132             :             }
     133             :         }
     134             : 
     135           1 :         if (GDALPamDataset::Close() != CE_None)
     136           0 :             eErr = CE_Failure;
     137             :     }
     138           2 :     return eErr;
     139             : }
     140             : 
     141             : /************************************************************************/
     142             : /*                          GetGeoTransform()                           */
     143             : /************************************************************************/
     144             : 
     145           1 : CPLErr DOQ2Dataset::GetGeoTransform(double *padfTransform)
     146             : 
     147             : {
     148           1 :     padfTransform[0] = dfULX;
     149           1 :     padfTransform[1] = dfXPixelSize;
     150           1 :     padfTransform[2] = 0.0;
     151           1 :     padfTransform[3] = dfULY;
     152           1 :     padfTransform[4] = 0.0;
     153           1 :     padfTransform[5] = -1 * dfYPixelSize;
     154             : 
     155           1 :     return CE_None;
     156             : }
     157             : 
     158             : /************************************************************************/
     159             : /*                                Open()                                */
     160             : /************************************************************************/
     161             : 
     162       29312 : GDALDataset *DOQ2Dataset::Open(GDALOpenInfo *poOpenInfo)
     163             : 
     164             : {
     165             :     /* -------------------------------------------------------------------- */
     166             :     /*      We assume the user is pointing to the binary (i.e. .bil) file.  */
     167             :     /* -------------------------------------------------------------------- */
     168       29312 :     if (poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fpL == nullptr)
     169       26516 :         return nullptr;
     170             : 
     171        2796 :     if (!STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
     172             :                         "BEGIN_USGS_DOQ_HEADER"))
     173        2795 :         return nullptr;
     174             : 
     175             :     /* -------------------------------------------------------------------- */
     176             :     /*      Confirm the requested access is supported.                      */
     177             :     /* -------------------------------------------------------------------- */
     178           1 :     if (poOpenInfo->eAccess == GA_Update)
     179             :     {
     180           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     181             :                  "The DOQ2 driver does not support update access to existing "
     182             :                  "datasets.");
     183           0 :         return nullptr;
     184             :     }
     185             : 
     186           1 :     int nBytesPerPixel = 0;
     187           1 :     int nWidth = 0;
     188           1 :     int nHeight = 0;
     189           1 :     int nBandStorage = 0;
     190           1 :     int nBandTypes = 0;
     191           1 :     const char *pszDatumLong = nullptr;
     192           1 :     const char *pszDatumShort = nullptr;
     193           1 :     const char *pszUnits = nullptr;
     194           1 :     int nZone = 0;
     195           1 :     int nProjType = 0;
     196           1 :     int nSkipBytes = 0;
     197           1 :     int nBandCount = 0;
     198           1 :     double dfULXMap = 0.0;
     199           1 :     double dfULYMap = 0.0;
     200           1 :     double dfXDim = 0.0;
     201           1 :     double dfYDim = 0.0;
     202           1 :     char **papszMetadata = nullptr;
     203             : 
     204             :     /* read and discard the first line */
     205           1 :     CPL_IGNORE_RET_VAL(CPLReadLineL(poOpenInfo->fpL));
     206             : 
     207           1 :     const char *pszLine = nullptr;
     208          46 :     while ((pszLine = CPLReadLineL(poOpenInfo->fpL)) != nullptr)
     209             :     {
     210          46 :         if (EQUAL(pszLine, "END_USGS_DOQ_HEADER"))
     211           0 :             break;
     212             : 
     213          46 :         char **papszTokens = CSLTokenizeString(pszLine);
     214          46 :         if (CSLCount(papszTokens) < 2)
     215             :         {
     216           1 :             CSLDestroy(papszTokens);
     217           1 :             break;
     218             :         }
     219             : 
     220          46 :         if (EQUAL(papszTokens[0], "SAMPLES_AND_LINES") &&
     221           1 :             CSLCount(papszTokens) >= 3)
     222             :         {
     223           1 :             nWidth = atoi(papszTokens[1]);
     224           1 :             nHeight = atoi(papszTokens[2]);
     225             :         }
     226          44 :         else if (EQUAL(papszTokens[0], "BYTE_COUNT"))
     227             :         {
     228           1 :             nSkipBytes = atoi(papszTokens[1]);
     229             :         }
     230          44 :         else if (EQUAL(papszTokens[0], "XY_ORIGIN") &&
     231           1 :                  CSLCount(papszTokens) >= 3)
     232             :         {
     233           1 :             dfULXMap = CPLAtof(papszTokens[1]);
     234           1 :             dfULYMap = CPLAtof(papszTokens[2]);
     235             :         }
     236          42 :         else if (EQUAL(papszTokens[0], "HORIZONTAL_RESOLUTION"))
     237             :         {
     238           1 :             dfXDim = CPLAtof(papszTokens[1]);
     239           1 :             dfYDim = dfXDim;
     240             :         }
     241          41 :         else if (EQUAL(papszTokens[0], "BAND_ORGANIZATION"))
     242             :         {
     243           1 :             if (EQUAL(papszTokens[1], "SINGLE FILE"))
     244           0 :                 nBandStorage = 1;
     245           1 :             if (EQUAL(papszTokens[1], "BSQ"))
     246           0 :                 nBandStorage = 1;
     247           1 :             if (EQUAL(papszTokens[1], "BIL"))
     248           0 :                 nBandStorage = 1;
     249           1 :             if (EQUAL(papszTokens[1], "BIP"))
     250           1 :                 nBandStorage = 4;
     251             :         }
     252          40 :         else if (EQUAL(papszTokens[0], "BAND_CONTENT"))
     253             :         {
     254           3 :             if (EQUAL(papszTokens[1], "BLACK&WHITE"))
     255           0 :                 nBandTypes = 1;
     256           3 :             else if (EQUAL(papszTokens[1], "COLOR"))
     257           0 :                 nBandTypes = 5;
     258           3 :             else if (EQUAL(papszTokens[1], "RGB"))
     259           0 :                 nBandTypes = 5;
     260           3 :             else if (EQUAL(papszTokens[1], "RED"))
     261           1 :                 nBandTypes = 5;
     262           2 :             else if (EQUAL(papszTokens[1], "GREEN"))
     263           1 :                 nBandTypes = 5;
     264           1 :             else if (EQUAL(papszTokens[1], "BLUE"))
     265           1 :                 nBandTypes = 5;
     266             : 
     267           3 :             nBandCount++;
     268             :         }
     269          37 :         else if (EQUAL(papszTokens[0], "BITS_PER_PIXEL"))
     270             :         {
     271           1 :             nBytesPerPixel = atoi(papszTokens[1]) / 8;
     272             :         }
     273          36 :         else if (EQUAL(papszTokens[0], "HORIZONTAL_COORDINATE_SYSTEM"))
     274             :         {
     275           1 :             if (EQUAL(papszTokens[1], "UTM"))
     276           1 :                 nProjType = 1;
     277           0 :             else if (EQUAL(papszTokens[1], "SPCS"))
     278           0 :                 nProjType = 2;
     279           0 :             else if (EQUAL(papszTokens[1], "GEOGRAPHIC"))
     280           0 :                 nProjType = 0;
     281             :         }
     282          35 :         else if (EQUAL(papszTokens[0], "COORDINATE_ZONE"))
     283             :         {
     284           1 :             nZone = atoi(papszTokens[1]);
     285             :         }
     286          34 :         else if (EQUAL(papszTokens[0], "HORIZONTAL_UNITS"))
     287             :         {
     288           1 :             if (EQUAL(papszTokens[1], "METERS"))
     289           1 :                 pszUnits = "UNIT[\"metre\",1]";
     290           0 :             else if (EQUAL(papszTokens[1], "FEET"))
     291           0 :                 pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
     292             :         }
     293          33 :         else if (EQUAL(papszTokens[0], "HORIZONTAL_DATUM"))
     294             :         {
     295           1 :             if (EQUAL(papszTokens[1], "NAD27"))
     296             :             {
     297           0 :                 pszDatumLong = NAD27_DATUM;
     298           0 :                 pszDatumShort = "NAD 27";
     299             :             }
     300           1 :             else if (EQUAL(papszTokens[1], " WGS72"))
     301             :             {
     302           0 :                 pszDatumLong = WGS72_DATUM;
     303           0 :                 pszDatumShort = "WGS 72";
     304             :             }
     305           1 :             else if (EQUAL(papszTokens[1], "WGS84"))
     306             :             {
     307           0 :                 pszDatumLong = WGS84_DATUM;
     308           0 :                 pszDatumShort = "WGS 84";
     309             :             }
     310           1 :             else if (EQUAL(papszTokens[1], "NAD83"))
     311             :             {
     312           1 :                 pszDatumLong = NAD83_DATUM;
     313           1 :                 pszDatumShort = "NAD 83";
     314             :             }
     315             :             else
     316             :             {
     317           0 :                 pszDatumLong = "DATUM[\"unknown\"]";
     318           0 :                 pszDatumShort = "unknown";
     319             :             }
     320             :         }
     321             :         else
     322             :         {
     323             :             /* we want to generically capture all the other metadata */
     324          32 :             CPLString osMetaDataValue;
     325             : 
     326         250 :             for (int iToken = 1; papszTokens[iToken] != nullptr; iToken++)
     327             :             {
     328         218 :                 if (EQUAL(papszTokens[iToken], "*"))
     329           1 :                     continue;
     330             : 
     331         217 :                 if (iToken > 1)
     332         186 :                     osMetaDataValue += " ";
     333         217 :                 osMetaDataValue += papszTokens[iToken];
     334             :             }
     335             :             papszMetadata =
     336          32 :                 CSLAddNameValue(papszMetadata, papszTokens[0], osMetaDataValue);
     337             :         }
     338             : 
     339          45 :         CSLDestroy(papszTokens);
     340             :     }
     341             : 
     342           1 :     CPLReadLineL(nullptr);
     343             : 
     344             :     /* -------------------------------------------------------------------- */
     345             :     /*      Do these values look coherent for a DOQ file?  It would be      */
     346             :     /*      nice to do a more comprehensive test than this!                 */
     347             :     /* -------------------------------------------------------------------- */
     348           1 :     if (nWidth < 500 || nWidth > 25000 || nHeight < 500 || nHeight > 25000 ||
     349           1 :         nBandStorage < 0 || nBandStorage > 4 || nBandTypes < 1 ||
     350           1 :         nBandTypes > 9 || nBytesPerPixel < 0)
     351             :     {
     352           0 :         CSLDestroy(papszMetadata);
     353           0 :         return nullptr;
     354             :     }
     355             : 
     356             :     /* -------------------------------------------------------------------- */
     357             :     /*      Check the configuration.  We don't currently handle all         */
     358             :     /*      variations, only the common ones.                               */
     359             :     /* -------------------------------------------------------------------- */
     360           1 :     if (nBandTypes > 5)
     361             :     {
     362           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     363             :                  "DOQ Data Type (%d) is not a supported configuration.",
     364             :                  nBandTypes);
     365           0 :         CSLDestroy(papszMetadata);
     366           0 :         return nullptr;
     367             :     }
     368             : 
     369             :     /* -------------------------------------------------------------------- */
     370             :     /*      Confirm the requested access is supported.                      */
     371             :     /* -------------------------------------------------------------------- */
     372           1 :     if (poOpenInfo->eAccess == GA_Update)
     373             :     {
     374           0 :         CSLDestroy(papszMetadata);
     375           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     376             :                  "The DOQ2 driver does not support update access to existing"
     377             :                  " datasets.");
     378           0 :         return nullptr;
     379             :     }
     380             :     /* -------------------------------------------------------------------- */
     381             :     /*      Create a corresponding GDALDataset.                             */
     382             :     /* -------------------------------------------------------------------- */
     383           2 :     auto poDS = std::make_unique<DOQ2Dataset>();
     384             : 
     385           1 :     poDS->nRasterXSize = nWidth;
     386           1 :     poDS->nRasterYSize = nHeight;
     387             : 
     388           1 :     poDS->SetMetadata(papszMetadata);
     389           1 :     CSLDestroy(papszMetadata);
     390           1 :     papszMetadata = nullptr;
     391             : 
     392           1 :     poDS->fpImage = poOpenInfo->fpL;
     393           1 :     poOpenInfo->fpL = nullptr;
     394             : 
     395             :     /* -------------------------------------------------------------------- */
     396             :     /*      Compute layout of data.                                         */
     397             :     /* -------------------------------------------------------------------- */
     398           1 :     if (nBandCount < 2)
     399             :     {
     400           0 :         nBandCount = nBytesPerPixel;
     401           0 :         if (!GDALCheckBandCount(nBandCount, FALSE))
     402             :         {
     403           0 :             return nullptr;
     404             :         }
     405             :     }
     406             :     else
     407             :     {
     408           1 :         if (nBytesPerPixel > INT_MAX / nBandCount)
     409             :         {
     410           0 :             return nullptr;
     411             :         }
     412           1 :         nBytesPerPixel *= nBandCount;
     413             :     }
     414             : 
     415           1 :     if (nBytesPerPixel > INT_MAX / nWidth)
     416             :     {
     417           0 :         return nullptr;
     418             :     }
     419           1 :     const int nBytesPerLine = nBytesPerPixel * nWidth;
     420             : 
     421             :     /* -------------------------------------------------------------------- */
     422             :     /*      Create band information objects.                                */
     423             :     /* -------------------------------------------------------------------- */
     424           4 :     for (int i = 0; i < nBandCount; i++)
     425             :     {
     426             :         auto poBand = RawRasterBand::Create(
     427           6 :             poDS.get(), i + 1, poDS->fpImage, nSkipBytes + i, nBytesPerPixel,
     428             :             nBytesPerLine, GDT_Byte,
     429             :             RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     430           3 :             RawRasterBand::OwnFP::NO);
     431           3 :         if (!poBand)
     432           0 :             return nullptr;
     433           3 :         poDS->SetBand(i + 1, std::move(poBand));
     434             :     }
     435             : 
     436           1 :     if (nProjType == 1)
     437             :     {
     438           2 :         poDS->m_oSRS.importFromWkt(
     439             :             CPLSPrintf(UTM_FORMAT, pszDatumShort ? pszDatumShort : "", nZone,
     440             :                        pszDatumLong ? pszDatumLong : "",
     441           1 :                        nZone >= 1 && nZone <= 60 ? nZone * 6 - 183 : 0,
     442             :                        pszUnits ? pszUnits : ""));
     443             :     }
     444             : 
     445           1 :     poDS->dfULX = dfULXMap;
     446           1 :     poDS->dfULY = dfULYMap;
     447             : 
     448           1 :     poDS->dfXPixelSize = dfXDim;
     449           1 :     poDS->dfYPixelSize = dfYDim;
     450             : 
     451             :     /* -------------------------------------------------------------------- */
     452             :     /*      Initialize any PAM information.                                 */
     453             :     /* -------------------------------------------------------------------- */
     454           1 :     poDS->SetDescription(poOpenInfo->pszFilename);
     455           1 :     poDS->TryLoadXML();
     456             : 
     457             :     /* -------------------------------------------------------------------- */
     458             :     /*      Check for overviews.                                            */
     459             :     /* -------------------------------------------------------------------- */
     460           1 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     461             : 
     462           1 :     return poDS.release();
     463             : }
     464             : 
     465             : /************************************************************************/
     466             : /*                         GDALRegister_DOQ1()                          */
     467             : /************************************************************************/
     468             : 
     469        1511 : void GDALRegister_DOQ2()
     470             : 
     471             : {
     472        1511 :     if (GDALGetDriverByName("DOQ2") != nullptr)
     473         295 :         return;
     474             : 
     475        1216 :     GDALDriver *poDriver = new GDALDriver();
     476             : 
     477        1216 :     poDriver->SetDescription("DOQ2");
     478        1216 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     479        1216 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "USGS DOQ (New Style)");
     480        1216 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/doq2.html");
     481        1216 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     482             : 
     483        1216 :     poDriver->pfnOpen = DOQ2Dataset::Open;
     484             : 
     485        1216 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     486             : }