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