Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  USGS DOQ Driver (First Generation Format)
       4             :  * Purpose:  Implementation of DOQ1Dataset
       5             :  * Author:   Frank Warmerdam, warmerda@home.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 1999, Frank Warmerdam
       9             :  * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "gdal_frmts.h"
      15             : #include "cpl_string.h"
      16             : #include "rawdataset.h"
      17             : 
      18             : #include <algorithm>
      19             : #include <cmath>
      20             : 
      21             : #ifndef UTM_FORMAT_defined
      22             : #define UTM_FORMAT_defined
      23             : 
      24             : static const char UTM_FORMAT[] =
      25             :     "PROJCS[\"%s / UTM zone %dN\",GEOGCS[%s,PRIMEM[\"Greenwich\",0],"
      26             :     "UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],"
      27             :     "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%d],"
      28             :     "PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],"
      29             :     "PARAMETER[\"false_northing\",0],%s]";
      30             : 
      31             : static const char WGS84_DATUM[] =
      32             :     "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]";
      33             : 
      34             : static const char WGS72_DATUM[] =
      35             :     "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]";
      36             : 
      37             : static const char NAD27_DATUM[] =
      38             :     "\"NAD27\",DATUM[\"North_American_Datum_1927\","
      39             :     "SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]";
      40             : 
      41             : static const char NAD83_DATUM[] =
      42             :     "\"NAD83\",DATUM[\"North_American_Datum_1983\","
      43             :     "SPHEROID[\"GRS 1980\",6378137,298.257222101]]";
      44             : 
      45             : #endif
      46             : 
      47             : /************************************************************************/
      48             : /*                            DOQGetField()                             */
      49             : /************************************************************************/
      50             : 
      51       11720 : static double DOQGetField(unsigned char *pabyData, int nBytes)
      52             : 
      53             : {
      54       11720 :     char szWork[128] = {'\0'};
      55             : 
      56       11720 :     memcpy(szWork, reinterpret_cast<const char *>(pabyData), nBytes);
      57       11720 :     szWork[nBytes] = '\0';
      58             : 
      59       64560 :     for (int i = 0; i < nBytes; i++)
      60             :     {
      61       52840 :         if (szWork[i] == 'D' || szWork[i] == 'd')
      62         433 :             szWork[i] = 'E';
      63             :     }
      64             : 
      65       23440 :     return CPLAtof(szWork);
      66             : }
      67             : 
      68             : /************************************************************************/
      69             : /*                         DOQGetDescription()                          */
      70             : /************************************************************************/
      71             : 
      72           2 : static void DOQGetDescription(GDALDataset *poDS, unsigned char *pabyData)
      73             : 
      74             : {
      75           2 :     char szWork[128] = {' '};
      76             : 
      77           2 :     const char *pszDescBegin = "USGS GeoTIFF DOQ 1:12000 Q-Quad of ";
      78           2 :     memcpy(szWork, pszDescBegin, strlen(pszDescBegin));
      79           2 :     memcpy(szWork + strlen(pszDescBegin),
      80             :            reinterpret_cast<const char *>(pabyData + 0), 38);
      81             : 
      82           2 :     int i = 0;
      83           2 :     while (*(szWork + 72 - i) == ' ')
      84             :     {
      85           0 :         i++;
      86             :     }
      87           2 :     i--;
      88             : 
      89           2 :     memcpy(szWork + 73 - i, reinterpret_cast<const char *>(pabyData + 38), 2);
      90           2 :     memcpy(szWork + 76 - i, reinterpret_cast<const char *>(pabyData + 44), 2);
      91           2 :     szWork[77 - i] = '\0';
      92             : 
      93           2 :     poDS->SetMetadataItem("DOQ_DESC", szWork);
      94           2 : }
      95             : 
      96             : /************************************************************************/
      97             : /* ==================================================================== */
      98             : /*                              DOQ1Dataset                             */
      99             : /* ==================================================================== */
     100             : /************************************************************************/
     101             : 
     102             : class DOQ1Dataset final : public RawDataset
     103             : {
     104             :     VSILFILE *fpImage = nullptr;  // Image data file.
     105             : 
     106             :     double dfULX = 0;
     107             :     double dfULY = 0;
     108             :     double dfXPixelSize = 0;
     109             :     double dfYPixelSize = 0;
     110             : 
     111             :     OGRSpatialReference m_oSRS{};
     112             : 
     113             :     CPL_DISALLOW_COPY_ASSIGN(DOQ1Dataset)
     114             : 
     115             :     CPLErr Close() override;
     116             : 
     117             :   public:
     118             :     DOQ1Dataset();
     119             :     ~DOQ1Dataset();
     120             : 
     121             :     CPLErr GetGeoTransform(double *padfTransform) override;
     122             : 
     123           0 :     const OGRSpatialReference *GetSpatialRef() const override
     124             :     {
     125           0 :         return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
     126             :     }
     127             : 
     128             :     static GDALDataset *Open(GDALOpenInfo *);
     129             : };
     130             : 
     131             : /************************************************************************/
     132             : /*                            DOQ1Dataset()                             */
     133             : /************************************************************************/
     134             : 
     135           2 : DOQ1Dataset::DOQ1Dataset()
     136             : {
     137           2 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     138           2 : }
     139             : 
     140             : /************************************************************************/
     141             : /*                            ~DOQ1Dataset()                            */
     142             : /************************************************************************/
     143             : 
     144           4 : DOQ1Dataset::~DOQ1Dataset()
     145             : 
     146             : {
     147           2 :     DOQ1Dataset::Close();
     148           4 : }
     149             : 
     150             : /************************************************************************/
     151             : /*                              Close()                                 */
     152             : /************************************************************************/
     153             : 
     154           4 : CPLErr DOQ1Dataset::Close()
     155             : {
     156           4 :     CPLErr eErr = CE_None;
     157           4 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     158             :     {
     159           2 :         if (DOQ1Dataset::FlushCache(true) != CE_None)
     160           0 :             eErr = CE_Failure;
     161             : 
     162           2 :         if (fpImage)
     163             :         {
     164           2 :             if (VSIFCloseL(fpImage) != 0)
     165             :             {
     166           0 :                 CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     167           0 :                 eErr = CE_Failure;
     168             :             }
     169             :         }
     170             : 
     171           2 :         if (GDALPamDataset::Close() != CE_None)
     172           0 :             eErr = CE_Failure;
     173             :     }
     174           4 :     return eErr;
     175             : }
     176             : 
     177             : /************************************************************************/
     178             : /*                          GetGeoTransform()                           */
     179             : /************************************************************************/
     180             : 
     181           0 : CPLErr DOQ1Dataset::GetGeoTransform(double *padfTransform)
     182             : 
     183             : {
     184           0 :     padfTransform[0] = dfULX;
     185           0 :     padfTransform[1] = dfXPixelSize;
     186           0 :     padfTransform[2] = 0.0;
     187           0 :     padfTransform[3] = dfULY;
     188           0 :     padfTransform[4] = 0.0;
     189           0 :     padfTransform[5] = -1 * dfYPixelSize;
     190             : 
     191           0 :     return CE_None;
     192             : }
     193             : 
     194             : /************************************************************************/
     195             : /*                                Open()                                */
     196             : /************************************************************************/
     197             : 
     198       31222 : GDALDataset *DOQ1Dataset::Open(GDALOpenInfo *poOpenInfo)
     199             : 
     200             : {
     201             :     /* -------------------------------------------------------------------- */
     202             :     /*      We assume the user is pointing to the binary (i.e. .bil) file.  */
     203             :     /* -------------------------------------------------------------------- */
     204       31222 :     if (poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fpL == nullptr)
     205       28296 :         return nullptr;
     206             : 
     207             :     /* -------------------------------------------------------------------- */
     208             :     /*      Attempt to extract a few key values from the header.            */
     209             :     /* -------------------------------------------------------------------- */
     210        2926 :     const double dfWidth = DOQGetField(poOpenInfo->pabyHeader + 150, 6);
     211        2926 :     const double dfHeight = DOQGetField(poOpenInfo->pabyHeader + 144, 6);
     212        2926 :     const double dfBandStorage = DOQGetField(poOpenInfo->pabyHeader + 162, 3);
     213        2926 :     const double dfBandTypes = DOQGetField(poOpenInfo->pabyHeader + 156, 3);
     214             : 
     215             :     /* -------------------------------------------------------------------- */
     216             :     /*      Do these values look coherent for a DOQ file?  It would be      */
     217             :     /*      nice to do a more comprehensive test than this!                 */
     218             :     /* -------------------------------------------------------------------- */
     219          50 :     if (dfWidth < 500 || dfWidth > 25000 || std::isnan(dfWidth) ||
     220           2 :         dfHeight < 500 || dfHeight > 25000 || std::isnan(dfHeight) ||
     221           2 :         dfBandStorage < 0 || dfBandStorage > 4 || std::isnan(dfBandStorage) ||
     222        2976 :         dfBandTypes < 1 || dfBandTypes > 9 || std::isnan(dfBandTypes))
     223        2924 :         return nullptr;
     224             : 
     225           2 :     const int nWidth = static_cast<int>(dfWidth);
     226           2 :     const int nHeight = static_cast<int>(dfHeight);
     227             :     /*const int nBandStorage = static_cast<int>(dfBandStorage);*/
     228           2 :     const int nBandTypes = static_cast<int>(dfBandTypes);
     229             : 
     230             :     /* -------------------------------------------------------------------- */
     231             :     /*      Check the configuration.  We don't currently handle all         */
     232             :     /*      variations, only the common ones.                               */
     233             :     /* -------------------------------------------------------------------- */
     234           2 :     if (nBandTypes > 5)
     235             :     {
     236           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     237             :                  "DOQ Data Type (%d) is not a supported configuration.",
     238             :                  nBandTypes);
     239           0 :         return nullptr;
     240             :     }
     241             : 
     242             :     /* -------------------------------------------------------------------- */
     243             :     /*      Confirm the requested access is supported.                      */
     244             :     /* -------------------------------------------------------------------- */
     245           2 :     if (poOpenInfo->eAccess == GA_Update)
     246             :     {
     247           0 :         ReportUpdateNotSupportedByDriver("DOQ1");
     248           0 :         return nullptr;
     249             :     }
     250             : 
     251             :     /* -------------------------------------------------------------------- */
     252             :     /*      Create a corresponding GDALDataset.                             */
     253             :     /* -------------------------------------------------------------------- */
     254           4 :     auto poDS = std::make_unique<DOQ1Dataset>();
     255             : 
     256             :     /* -------------------------------------------------------------------- */
     257             :     /*      Capture some information from the file that is of interest.     */
     258             :     /* -------------------------------------------------------------------- */
     259           2 :     poDS->nRasterXSize = nWidth;
     260           2 :     poDS->nRasterYSize = nHeight;
     261             : 
     262           2 :     std::swap(poDS->fpImage, poOpenInfo->fpL);
     263             : 
     264             :     /* -------------------------------------------------------------------- */
     265             :     /*      Compute layout of data.                                         */
     266             :     /* -------------------------------------------------------------------- */
     267           2 :     int nBytesPerPixel = 0;
     268             : 
     269           2 :     if (nBandTypes < 5)
     270           2 :         nBytesPerPixel = 1;
     271             :     else /* if( nBandTypes == 5 ) */
     272           0 :         nBytesPerPixel = 3;
     273             : 
     274           2 :     const int nBytesPerLine = nBytesPerPixel * nWidth;
     275           2 :     const int nSkipBytes = 4 * nBytesPerLine;
     276             : 
     277             :     /* -------------------------------------------------------------------- */
     278             :     /*      Create band information objects.                                */
     279             :     /* -------------------------------------------------------------------- */
     280           4 :     for (int i = 0; i < nBytesPerPixel; i++)
     281             :     {
     282             :         auto poBand = RawRasterBand::Create(
     283           4 :             poDS.get(), i + 1, poDS->fpImage, nSkipBytes + i, nBytesPerPixel,
     284             :             nBytesPerLine, GDT_Byte,
     285             :             RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     286           2 :             RawRasterBand::OwnFP::NO);
     287           2 :         if (!poBand)
     288           0 :             return nullptr;
     289           2 :         poDS->SetBand(i + 1, std::move(poBand));
     290             :     }
     291             : 
     292             :     /* -------------------------------------------------------------------- */
     293             :     /*      Set the description.                                            */
     294             :     /* -------------------------------------------------------------------- */
     295           2 :     DOQGetDescription(poDS.get(), poOpenInfo->pabyHeader);
     296             : 
     297             :     /* -------------------------------------------------------------------- */
     298             :     /*      Establish the projection string.                                */
     299             :     /* -------------------------------------------------------------------- */
     300           2 :     if (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 195, 3)) == 1)
     301             :     {
     302             :         int nZone =
     303           2 :             static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 198, 6));
     304           2 :         if (nZone < 0 || nZone > 60)
     305           0 :             nZone = 0;
     306             : 
     307           2 :         const char *pszUnits = nullptr;
     308           2 :         if (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 204, 3)) == 1)
     309           0 :             pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
     310             :         else
     311           2 :             pszUnits = "UNIT[\"metre\",1]";
     312             : 
     313           2 :         const char *pszDatumLong = nullptr;
     314           2 :         const char *pszDatumShort = nullptr;
     315           2 :         switch (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 167, 2)))
     316             :         {
     317           0 :             case 1:
     318           0 :                 pszDatumLong = NAD27_DATUM;
     319           0 :                 pszDatumShort = "NAD 27";
     320           0 :                 break;
     321             : 
     322           0 :             case 2:
     323           0 :                 pszDatumLong = WGS72_DATUM;
     324           0 :                 pszDatumShort = "WGS 72";
     325           0 :                 break;
     326             : 
     327           2 :             case 3:
     328           2 :                 pszDatumLong = WGS84_DATUM;
     329           2 :                 pszDatumShort = "WGS 84";
     330           2 :                 break;
     331             : 
     332           0 :             case 4:
     333           0 :                 pszDatumLong = NAD83_DATUM;
     334           0 :                 pszDatumShort = "NAD 83";
     335           0 :                 break;
     336             : 
     337           0 :             default:
     338           0 :                 pszDatumLong = "DATUM[\"unknown\"]";
     339           0 :                 pszDatumShort = "unknown";
     340           0 :                 break;
     341             :         }
     342             : 
     343           4 :         poDS->m_oSRS.importFromWkt(CPLSPrintf(UTM_FORMAT, pszDatumShort, nZone,
     344           2 :                                               pszDatumLong, nZone * 6 - 183,
     345             :                                               pszUnits));
     346             :     }
     347             : 
     348             :     /* -------------------------------------------------------------------- */
     349             :     /*      Read the georeferencing information.                            */
     350             :     /* -------------------------------------------------------------------- */
     351           2 :     unsigned char abyRecordData[500] = {'\0'};
     352             : 
     353           4 :     if (VSIFSeekL(poDS->fpImage, nBytesPerLine * 2, SEEK_SET) != 0 ||
     354           2 :         VSIFReadL(abyRecordData, sizeof(abyRecordData), 1, poDS->fpImage) != 1)
     355             :     {
     356           0 :         CPLError(CE_Failure, CPLE_FileIO, "Header read error on %s.",
     357             :                  poOpenInfo->pszFilename);
     358           0 :         return nullptr;
     359             :     }
     360             : 
     361           2 :     poDS->dfULX = DOQGetField(abyRecordData + 288, 24);
     362           2 :     poDS->dfULY = DOQGetField(abyRecordData + 312, 24);
     363             : 
     364           4 :     if (VSIFSeekL(poDS->fpImage, nBytesPerLine * 3, SEEK_SET) != 0 ||
     365           2 :         VSIFReadL(abyRecordData, sizeof(abyRecordData), 1, poDS->fpImage) != 1)
     366             :     {
     367           0 :         CPLError(CE_Failure, CPLE_FileIO, "Header read error on %s.",
     368             :                  poOpenInfo->pszFilename);
     369           0 :         return nullptr;
     370             :     }
     371             : 
     372           2 :     poDS->dfXPixelSize = DOQGetField(abyRecordData + 59, 12);
     373           2 :     poDS->dfYPixelSize = DOQGetField(abyRecordData + 71, 12);
     374             : 
     375             :     /* -------------------------------------------------------------------- */
     376             :     /*      Initialize any PAM information.                                 */
     377             :     /* -------------------------------------------------------------------- */
     378           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     379           2 :     poDS->TryLoadXML();
     380             : 
     381             :     /* -------------------------------------------------------------------- */
     382             :     /*      Check for overviews.                                            */
     383             :     /* -------------------------------------------------------------------- */
     384           2 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     385             : 
     386           2 :     return poDS.release();
     387             : }
     388             : 
     389             : /************************************************************************/
     390             : /*                         GDALRegister_DOQ1()                          */
     391             : /************************************************************************/
     392             : 
     393        1686 : void GDALRegister_DOQ1()
     394             : 
     395             : {
     396        1686 :     if (GDALGetDriverByName("DOQ1") != nullptr)
     397         302 :         return;
     398             : 
     399        1384 :     GDALDriver *poDriver = new GDALDriver();
     400             : 
     401        1384 :     poDriver->SetDescription("DOQ1");
     402        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     403        1384 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "USGS DOQ (Old Style)");
     404        1384 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/doq1.html");
     405        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     406             : 
     407        1384 :     poDriver->pfnOpen = DOQ1Dataset::Open;
     408             : 
     409        1384 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     410             : }