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       11512 : static double DOQGetField(unsigned char *pabyData, int nBytes)
      52             : 
      53             : {
      54       11512 :     char szWork[128] = {'\0'};
      55             : 
      56       11512 :     memcpy(szWork, reinterpret_cast<const char *>(pabyData), nBytes);
      57       11512 :     szWork[nBytes] = '\0';
      58             : 
      59       63416 :     for (int i = 0; i < nBytes; i++)
      60             :     {
      61       51904 :         if (szWork[i] == 'D' || szWork[i] == 'd')
      62         436 :             szWork[i] = 'E';
      63             :     }
      64             : 
      65       23024 :     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       30750 : GDALDataset *DOQ1Dataset::Open(GDALOpenInfo *poOpenInfo)
     199             : 
     200             : {
     201             :     /* -------------------------------------------------------------------- */
     202             :     /*      We assume the user is pointing to the binary (i.e. .bil) file.  */
     203             :     /* -------------------------------------------------------------------- */
     204       30750 :     if (poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fpL == nullptr)
     205       27876 :         return nullptr;
     206             : 
     207             :     /* -------------------------------------------------------------------- */
     208             :     /*      Attempt to extract a few key values from the header.            */
     209             :     /* -------------------------------------------------------------------- */
     210        2874 :     const double dfWidth = DOQGetField(poOpenInfo->pabyHeader + 150, 6);
     211        2874 :     const double dfHeight = DOQGetField(poOpenInfo->pabyHeader + 144, 6);
     212        2874 :     const double dfBandStorage = DOQGetField(poOpenInfo->pabyHeader + 162, 3);
     213        2874 :     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          52 :     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        2926 :         dfBandTypes < 1 || dfBandTypes > 9 || std::isnan(dfBandTypes))
     223        2872 :         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 :         CPLError(CE_Failure, CPLE_NotSupported,
     248             :                  "The DOQ1 driver does not support update access to existing "
     249             :                  "datasets.");
     250           0 :         return nullptr;
     251             :     }
     252             : 
     253             :     /* -------------------------------------------------------------------- */
     254             :     /*      Create a corresponding GDALDataset.                             */
     255             :     /* -------------------------------------------------------------------- */
     256           4 :     auto poDS = std::make_unique<DOQ1Dataset>();
     257             : 
     258             :     /* -------------------------------------------------------------------- */
     259             :     /*      Capture some information from the file that is of interest.     */
     260             :     /* -------------------------------------------------------------------- */
     261           2 :     poDS->nRasterXSize = nWidth;
     262           2 :     poDS->nRasterYSize = nHeight;
     263             : 
     264           2 :     std::swap(poDS->fpImage, poOpenInfo->fpL);
     265             : 
     266             :     /* -------------------------------------------------------------------- */
     267             :     /*      Compute layout of data.                                         */
     268             :     /* -------------------------------------------------------------------- */
     269           2 :     int nBytesPerPixel = 0;
     270             : 
     271           2 :     if (nBandTypes < 5)
     272           2 :         nBytesPerPixel = 1;
     273             :     else /* if( nBandTypes == 5 ) */
     274           0 :         nBytesPerPixel = 3;
     275             : 
     276           2 :     const int nBytesPerLine = nBytesPerPixel * nWidth;
     277           2 :     const int nSkipBytes = 4 * nBytesPerLine;
     278             : 
     279             :     /* -------------------------------------------------------------------- */
     280             :     /*      Create band information objects.                                */
     281             :     /* -------------------------------------------------------------------- */
     282           4 :     for (int i = 0; i < nBytesPerPixel; i++)
     283             :     {
     284             :         auto poBand = RawRasterBand::Create(
     285           4 :             poDS.get(), i + 1, poDS->fpImage, nSkipBytes + i, nBytesPerPixel,
     286             :             nBytesPerLine, GDT_Byte,
     287             :             RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     288           2 :             RawRasterBand::OwnFP::NO);
     289           2 :         if (!poBand)
     290           0 :             return nullptr;
     291           2 :         poDS->SetBand(i + 1, std::move(poBand));
     292             :     }
     293             : 
     294             :     /* -------------------------------------------------------------------- */
     295             :     /*      Set the description.                                            */
     296             :     /* -------------------------------------------------------------------- */
     297           2 :     DOQGetDescription(poDS.get(), poOpenInfo->pabyHeader);
     298             : 
     299             :     /* -------------------------------------------------------------------- */
     300             :     /*      Establish the projection string.                                */
     301             :     /* -------------------------------------------------------------------- */
     302           2 :     if (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 195, 3)) == 1)
     303             :     {
     304             :         int nZone =
     305           2 :             static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 198, 6));
     306           2 :         if (nZone < 0 || nZone > 60)
     307           0 :             nZone = 0;
     308             : 
     309           2 :         const char *pszUnits = nullptr;
     310           2 :         if (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 204, 3)) == 1)
     311           0 :             pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
     312             :         else
     313           2 :             pszUnits = "UNIT[\"metre\",1]";
     314             : 
     315           2 :         const char *pszDatumLong = nullptr;
     316           2 :         const char *pszDatumShort = nullptr;
     317           2 :         switch (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 167, 2)))
     318             :         {
     319           0 :             case 1:
     320           0 :                 pszDatumLong = NAD27_DATUM;
     321           0 :                 pszDatumShort = "NAD 27";
     322           0 :                 break;
     323             : 
     324           0 :             case 2:
     325           0 :                 pszDatumLong = WGS72_DATUM;
     326           0 :                 pszDatumShort = "WGS 72";
     327           0 :                 break;
     328             : 
     329           2 :             case 3:
     330           2 :                 pszDatumLong = WGS84_DATUM;
     331           2 :                 pszDatumShort = "WGS 84";
     332           2 :                 break;
     333             : 
     334           0 :             case 4:
     335           0 :                 pszDatumLong = NAD83_DATUM;
     336           0 :                 pszDatumShort = "NAD 83";
     337           0 :                 break;
     338             : 
     339           0 :             default:
     340           0 :                 pszDatumLong = "DATUM[\"unknown\"]";
     341           0 :                 pszDatumShort = "unknown";
     342           0 :                 break;
     343             :         }
     344             : 
     345           4 :         poDS->m_oSRS.importFromWkt(CPLSPrintf(UTM_FORMAT, pszDatumShort, nZone,
     346           2 :                                               pszDatumLong, nZone * 6 - 183,
     347             :                                               pszUnits));
     348             :     }
     349             : 
     350             :     /* -------------------------------------------------------------------- */
     351             :     /*      Read the georeferencing information.                            */
     352             :     /* -------------------------------------------------------------------- */
     353           2 :     unsigned char abyRecordData[500] = {'\0'};
     354             : 
     355           4 :     if (VSIFSeekL(poDS->fpImage, nBytesPerLine * 2, SEEK_SET) != 0 ||
     356           2 :         VSIFReadL(abyRecordData, sizeof(abyRecordData), 1, poDS->fpImage) != 1)
     357             :     {
     358           0 :         CPLError(CE_Failure, CPLE_FileIO, "Header read error on %s.",
     359             :                  poOpenInfo->pszFilename);
     360           0 :         return nullptr;
     361             :     }
     362             : 
     363           2 :     poDS->dfULX = DOQGetField(abyRecordData + 288, 24);
     364           2 :     poDS->dfULY = DOQGetField(abyRecordData + 312, 24);
     365             : 
     366           4 :     if (VSIFSeekL(poDS->fpImage, nBytesPerLine * 3, SEEK_SET) != 0 ||
     367           2 :         VSIFReadL(abyRecordData, sizeof(abyRecordData), 1, poDS->fpImage) != 1)
     368             :     {
     369           0 :         CPLError(CE_Failure, CPLE_FileIO, "Header read error on %s.",
     370             :                  poOpenInfo->pszFilename);
     371           0 :         return nullptr;
     372             :     }
     373             : 
     374           2 :     poDS->dfXPixelSize = DOQGetField(abyRecordData + 59, 12);
     375           2 :     poDS->dfYPixelSize = DOQGetField(abyRecordData + 71, 12);
     376             : 
     377             :     /* -------------------------------------------------------------------- */
     378             :     /*      Initialize any PAM information.                                 */
     379             :     /* -------------------------------------------------------------------- */
     380           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     381           2 :     poDS->TryLoadXML();
     382             : 
     383             :     /* -------------------------------------------------------------------- */
     384             :     /*      Check for overviews.                                            */
     385             :     /* -------------------------------------------------------------------- */
     386           2 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     387             : 
     388           2 :     return poDS.release();
     389             : }
     390             : 
     391             : /************************************************************************/
     392             : /*                         GDALRegister_DOQ1()                          */
     393             : /************************************************************************/
     394             : 
     395        1595 : void GDALRegister_DOQ1()
     396             : 
     397             : {
     398        1595 :     if (GDALGetDriverByName("DOQ1") != nullptr)
     399         302 :         return;
     400             : 
     401        1293 :     GDALDriver *poDriver = new GDALDriver();
     402             : 
     403        1293 :     poDriver->SetDescription("DOQ1");
     404        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     405        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "USGS DOQ (Old Style)");
     406        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/doq1.html");
     407        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     408             : 
     409        1293 :     poDriver->pfnOpen = DOQ1Dataset::Open;
     410             : 
     411        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     412             : }