Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  R Format Driver
       4             :  * Purpose:  Read/write R stats package object format.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2009, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2009-2010, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "rdataset.h"
      16             : 
      17             : #include <cstddef>
      18             : #include <cstdlib>
      19             : #include <cstring>
      20             : #if HAVE_FCNTL_H
      21             : #include <fcntl.h>
      22             : #endif
      23             : 
      24             : #include <algorithm>
      25             : #include <limits>
      26             : #include <string>
      27             : #include <utility>
      28             : 
      29             : #include "cpl_conv.h"
      30             : #include "cpl_error.h"
      31             : #include "cpl_progress.h"
      32             : #include "cpl_string.h"
      33             : #include "cpl_vsi.h"
      34             : #include "gdal.h"
      35             : #include "gdal_frmts.h"
      36             : #include "gdal_pam.h"
      37             : #include "gdal_priv.h"
      38             : 
      39             : // constexpr int R_NILSXP = 0;
      40             : constexpr int R_LISTSXP = 2;
      41             : constexpr int R_CHARSXP = 9;
      42             : constexpr int R_INTSXP = 13;
      43             : constexpr int R_REALSXP = 14;
      44             : constexpr int R_STRSXP = 16;
      45             : 
      46             : namespace
      47             : {
      48             : 
      49             : // TODO(schwehr): Move this to port/? for general use.
      50          20 : bool SafeMult(GIntBig a, GIntBig b, GIntBig *result)
      51             : {
      52          20 :     if (a == 0 || b == 0)
      53             :     {
      54           0 :         *result = 0;
      55           0 :         return true;
      56             :     }
      57             : 
      58          20 :     bool result_positive = (a >= 0 && b >= 0) || (a < 0 && b < 0);
      59          20 :     if (result_positive)
      60             :     {
      61             :         // Cannot convert min() to positive.
      62          40 :         if (a == std::numeric_limits<GIntBig>::min() ||
      63          20 :             b == std::numeric_limits<GIntBig>::min())
      64             :         {
      65           0 :             *result = 0;
      66           0 :             return false;
      67             :         }
      68          20 :         if (a < 0)
      69             :         {
      70           0 :             a = -a;
      71           0 :             b = -b;
      72             :         }
      73          20 :         if (a > std::numeric_limits<GIntBig>::max() / b)
      74             :         {
      75           0 :             *result = 0;
      76           0 :             return false;
      77             :         }
      78          20 :         *result = a * b;
      79          20 :         return true;
      80             :     }
      81             : 
      82           0 :     if (b < a)
      83           0 :         std::swap(a, b);
      84           0 :     if (a < (std::numeric_limits<GIntBig>::min() + 1) / b)
      85             :     {
      86           0 :         *result = 0;
      87           0 :         return false;
      88             :     }
      89             : 
      90           0 :     *result = a * b;
      91           0 :     return true;
      92             : }
      93             : 
      94             : }  // namespace
      95             : 
      96             : /************************************************************************/
      97             : /*                            RRasterBand()                             */
      98             : /************************************************************************/
      99             : 
     100           7 : RRasterBand::RRasterBand(RDataset *poDSIn, int nBandIn,
     101           7 :                          const double *padfMatrixValuesIn)
     102           7 :     : padfMatrixValues(padfMatrixValuesIn)
     103             : {
     104           7 :     poDS = poDSIn;
     105           7 :     nBand = nBandIn;
     106             : 
     107           7 :     eDataType = GDT_Float64;
     108             : 
     109           7 :     nBlockXSize = poDSIn->nRasterXSize;
     110           7 :     nBlockYSize = 1;
     111           7 : }
     112             : 
     113             : /************************************************************************/
     114             : /*                             IReadBlock()                             */
     115             : /************************************************************************/
     116             : 
     117          45 : CPLErr RRasterBand::IReadBlock(int /* nBlockXOff */, int nBlockYOff,
     118             :                                void *pImage)
     119             : {
     120          45 :     memcpy(pImage, padfMatrixValues + nBlockYOff * nBlockXSize,
     121          45 :            nBlockXSize * 8);
     122          45 :     return CE_None;
     123             : }
     124             : 
     125             : /************************************************************************/
     126             : /* ==================================================================== */
     127             : /*                              RDataset()                              */
     128             : /* ==================================================================== */
     129             : /************************************************************************/
     130             : 
     131             : /************************************************************************/
     132             : /*                              RDataset()                              */
     133             : /************************************************************************/
     134             : 
     135          10 : RDataset::RDataset()
     136          10 :     : fp(nullptr), bASCII(FALSE), nStartOfData(0), padfMatrixValues(nullptr)
     137             : {
     138          10 : }
     139             : 
     140             : /************************************************************************/
     141             : /*                             ~RDataset()                              */
     142             : /************************************************************************/
     143             : 
     144          20 : RDataset::~RDataset()
     145             : {
     146          10 :     FlushCache(true);
     147          10 :     CPLFree(padfMatrixValues);
     148             : 
     149          10 :     if (fp)
     150          10 :         VSIFCloseL(fp);
     151          20 : }
     152             : 
     153             : /************************************************************************/
     154             : /*                             ASCIIFGets()                             */
     155             : /*                                                                      */
     156             : /*      Fetch one line from an ASCII source into osLastStringRead.      */
     157             : /************************************************************************/
     158             : 
     159        1385 : const char *RDataset::ASCIIFGets()
     160             : 
     161             : {
     162        1385 :     char chNextChar = '\0';
     163             : 
     164        1385 :     osLastStringRead.resize(0);
     165             : 
     166        3952 :     do
     167             :     {
     168        5337 :         chNextChar = '\n';
     169        5337 :         VSIFReadL(&chNextChar, 1, 1, fp);
     170        5337 :         if (chNextChar != '\n')
     171        3952 :             osLastStringRead += chNextChar;
     172        5337 :     } while (chNextChar != '\n' && chNextChar != '\0');
     173             : 
     174        2770 :     return osLastStringRead;
     175             : }
     176             : 
     177             : /************************************************************************/
     178             : /*                            ReadInteger()                             */
     179             : /************************************************************************/
     180             : 
     181         190 : int RDataset::ReadInteger()
     182             : 
     183             : {
     184         190 :     if (bASCII)
     185             :     {
     186          95 :         return atoi(ASCIIFGets());
     187             :     }
     188             : 
     189          95 :     GInt32 nValue = 0;
     190             : 
     191          95 :     if (VSIFReadL(&nValue, 4, 1, fp) != 1)
     192           0 :         return -1;
     193          95 :     CPL_MSBPTR32(&nValue);
     194             : 
     195          95 :     return nValue;
     196             : }
     197             : 
     198             : /************************************************************************/
     199             : /*                             ReadFloat()                              */
     200             : /************************************************************************/
     201             : 
     202        1280 : double RDataset::ReadFloat()
     203             : 
     204             : {
     205        1280 :     if (bASCII)
     206             :     {
     207        1280 :         return CPLAtof(ASCIIFGets());
     208             :     }
     209             : 
     210           0 :     double dfValue = 0.0;
     211             : 
     212           0 :     if (VSIFReadL(&dfValue, 8, 1, fp) != 1)
     213           0 :         return -1;
     214           0 :     CPL_MSBPTR64(&dfValue);
     215             : 
     216           0 :     return dfValue;
     217             : }
     218             : 
     219             : /************************************************************************/
     220             : /*                             ReadString()                             */
     221             : /************************************************************************/
     222             : 
     223          20 : const char *RDataset::ReadString()
     224             : 
     225             : {
     226          20 :     if (ReadInteger() % 256 != R_CHARSXP)
     227             :     {
     228           0 :         osLastStringRead = "";
     229           0 :         return "";
     230             :     }
     231             : 
     232          20 :     const int nLenSigned = ReadInteger();
     233          20 :     if (nLenSigned < 0)
     234             :     {
     235           0 :         osLastStringRead = "";
     236           0 :         return "";
     237             :     }
     238          20 :     const size_t nLen = static_cast<size_t>(nLenSigned);
     239             : 
     240          20 :     char *pachWrkBuf = static_cast<char *>(VSIMalloc(nLen));
     241          20 :     if (pachWrkBuf == nullptr)
     242             :     {
     243           0 :         osLastStringRead = "";
     244           0 :         return "";
     245             :     }
     246          20 :     if (VSIFReadL(pachWrkBuf, 1, nLen, fp) != nLen)
     247             :     {
     248           0 :         osLastStringRead = "";
     249           0 :         CPLFree(pachWrkBuf);
     250           0 :         return "";
     251             :     }
     252             : 
     253          20 :     if (bASCII)
     254             :     {
     255             :         // Suck up newline and any extra junk.
     256          10 :         ASCIIFGets();
     257             :     }
     258             : 
     259          20 :     osLastStringRead.assign(pachWrkBuf, nLen);
     260          20 :     CPLFree(pachWrkBuf);
     261             : 
     262          20 :     return osLastStringRead;
     263             : }
     264             : 
     265             : /************************************************************************/
     266             : /*                              ReadPair()                              */
     267             : /************************************************************************/
     268             : 
     269          30 : bool RDataset::ReadPair(CPLString &osObjName, int &nObjCode)
     270             : 
     271             : {
     272          30 :     nObjCode = ReadInteger();
     273          30 :     if (nObjCode == 254)
     274          10 :         return true;
     275             : 
     276          20 :     if ((nObjCode % 256) != R_LISTSXP)
     277             :     {
     278           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     279             :                  "Did not find expected object pair object.");
     280           0 :         return false;
     281             :     }
     282             : 
     283          20 :     int nPairCount = ReadInteger();
     284          20 :     if (nPairCount != 1)
     285             :     {
     286           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     287             :                  "Did not find expected pair count of 1.");
     288           0 :         return false;
     289             :     }
     290             : 
     291             :     // Read the object name.
     292          20 :     const char *pszName = ReadString();
     293          20 :     if (pszName == nullptr || pszName[0] == '\0')
     294           0 :         return false;
     295             : 
     296          20 :     osObjName = pszName;
     297             : 
     298             :     // Confirm that we have a numeric matrix object.
     299          20 :     nObjCode = ReadInteger();
     300             : 
     301          20 :     return true;
     302             : }
     303             : 
     304             : /************************************************************************/
     305             : /*                              Identify()                              */
     306             : /************************************************************************/
     307             : 
     308       53821 : int RDataset::Identify(GDALOpenInfo *poOpenInfo)
     309             : {
     310       53821 :     if (poOpenInfo->nHeaderBytes < 50)
     311       49502 :         return FALSE;
     312             : 
     313             :     // If the extension is .rda and the file type is gzip
     314             :     // compressed we assume it is a gzipped R binary file.
     315        4353 :     if (memcmp(poOpenInfo->pabyHeader, "\037\213\b", 3) == 0 &&
     316          34 :         poOpenInfo->IsExtensionEqualToCI("rda"))
     317           6 :         return TRUE;
     318             : 
     319             :     // Is this an ASCII or XDR binary R file?
     320        4313 :     if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "RDA2\nA\n") &&
     321        4303 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "RDX2\nX\n"))
     322        4299 :         return FALSE;
     323             : 
     324          14 :     return TRUE;
     325             : }
     326             : 
     327             : /************************************************************************/
     328             : /*                                Open()                                */
     329             : /************************************************************************/
     330             : 
     331          10 : GDALDataset *RDataset::Open(GDALOpenInfo *poOpenInfo)
     332             : {
     333             : #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
     334             :     if (poOpenInfo->pabyHeader == nullptr)
     335             :         return nullptr;
     336             : #else
     337             :     // During fuzzing, do not use Identify to reject crazy content.
     338          10 :     if (!Identify(poOpenInfo))
     339           0 :         return nullptr;
     340             : #endif
     341             : 
     342             :     // Confirm the requested access is supported.
     343          10 :     if (poOpenInfo->eAccess == GA_Update)
     344             :     {
     345           0 :         ReportUpdateNotSupportedByDriver("R");
     346           0 :         return nullptr;
     347             :     }
     348             : 
     349             :     // Do we need to route the file through the decompression machinery?
     350          10 :     const bool bCompressed =
     351          10 :         memcmp(poOpenInfo->pabyHeader, "\037\213\b", 3) == 0;
     352             :     const CPLString osAdjustedFilename =
     353          30 :         std::string(bCompressed ? "/vsigzip/" : "") + poOpenInfo->pszFilename;
     354             : 
     355             :     // Establish this as a dataset and open the file using VSI*L.
     356          20 :     auto poDS = std::make_unique<RDataset>();
     357             : 
     358          10 :     poDS->fp = VSIFOpenL(osAdjustedFilename, "r");
     359          10 :     if (poDS->fp == nullptr)
     360             :     {
     361           0 :         return nullptr;
     362             :     }
     363             : 
     364          10 :     poDS->bASCII = STARTS_WITH_CI(
     365             :         reinterpret_cast<char *>(poOpenInfo->pabyHeader), "RDA2\nA\n");
     366             : 
     367             :     // Confirm this is a version 2 file.
     368          10 :     VSIFSeekL(poDS->fp, 7, SEEK_SET);
     369          10 :     if (poDS->ReadInteger() != R_LISTSXP)
     370             :     {
     371           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     372             :                  "It appears %s is not a version 2 R object file after all!",
     373             :                  poOpenInfo->pszFilename);
     374           0 :         return nullptr;
     375             :     }
     376             : 
     377             :     // Skip the version values.
     378          10 :     poDS->ReadInteger();
     379          10 :     poDS->ReadInteger();
     380             : 
     381             :     // Confirm we have a numeric vector object in a pairlist.
     382          20 :     CPLString osObjName;
     383          10 :     int nObjCode = 0;
     384             : 
     385          10 :     if (!poDS->ReadPair(osObjName, nObjCode))
     386             :     {
     387           0 :         return nullptr;
     388             :     }
     389             : 
     390          10 :     if (nObjCode % 256 != R_REALSXP)
     391             :     {
     392           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     393             :                  "Failed to find expected numeric vector object.");
     394           0 :         return nullptr;
     395             :     }
     396             : 
     397          10 :     poDS->SetMetadataItem("R_OBJECT_NAME", osObjName);
     398             : 
     399             :     // Read the count.
     400          10 :     const int nValueCount = poDS->ReadInteger();
     401          10 :     if (nValueCount < 0)
     402             :     {
     403           0 :         CPLError(CE_Failure, CPLE_AppDefined, "nValueCount < 0: %d",
     404             :                  nValueCount);
     405           0 :         return nullptr;
     406             :     }
     407             : 
     408          10 :     poDS->nStartOfData = VSIFTellL(poDS->fp);
     409             : 
     410             :     // TODO(schwehr): Factor in the size of doubles.
     411             :     VSIStatBufL stat;
     412             :     const int dStatSuccess =
     413          10 :         VSIStatExL(osAdjustedFilename, &stat, VSI_STAT_SIZE_FLAG);
     414          20 :     if (dStatSuccess != 0 || static_cast<vsi_l_offset>(nValueCount) >
     415          10 :                                  stat.st_size - poDS->nStartOfData)
     416             :     {
     417           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     418             :                  "Corrupt file.  "
     419             :                  "Object claims to be larger than available bytes. "
     420             :                  "%d > " CPL_FRMT_GUIB,
     421           0 :                  nValueCount, stat.st_size - poDS->nStartOfData);
     422           0 :         return nullptr;
     423             :     }
     424             : 
     425             :     // Read/Skip ahead to attributes.
     426          10 :     if (poDS->bASCII)
     427             :     {
     428          10 :         poDS->padfMatrixValues =
     429           5 :             static_cast<double *>(VSIMalloc2(nValueCount, sizeof(double)));
     430           5 :         if (poDS->padfMatrixValues == nullptr)
     431             :         {
     432           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot allocate %d doubles",
     433             :                      nValueCount);
     434           0 :             return nullptr;
     435             :         }
     436        1285 :         for (int iValue = 0; iValue < nValueCount; iValue++)
     437        1280 :             poDS->padfMatrixValues[iValue] = poDS->ReadFloat();
     438             :     }
     439             :     else
     440             :     {
     441           5 :         VSIFSeekL(poDS->fp, 8 * nValueCount, SEEK_CUR);
     442             :     }
     443             : 
     444             :     // Read pairs till we run out, trying to find a few items that
     445             :     // have special meaning to us.
     446          10 :     poDS->nRasterXSize = 0;
     447          10 :     poDS->nRasterYSize = 0;
     448          10 :     int nBandCount = 0;
     449             : 
     450          20 :     while (poDS->ReadPair(osObjName, nObjCode) && nObjCode != 254)
     451             :     {
     452          10 :         if (osObjName == "dim" && nObjCode % 256 == R_INTSXP)
     453             :         {
     454          10 :             const int nCount = poDS->ReadInteger();
     455          10 :             if (nCount == 2)
     456             :             {
     457           0 :                 poDS->nRasterXSize = poDS->ReadInteger();
     458           0 :                 poDS->nRasterYSize = poDS->ReadInteger();
     459           0 :                 nBandCount = 1;
     460             :             }
     461          10 :             else if (nCount == 3)
     462             :             {
     463          10 :                 poDS->nRasterXSize = poDS->ReadInteger();
     464          10 :                 poDS->nRasterYSize = poDS->ReadInteger();
     465          10 :                 nBandCount = poDS->ReadInteger();
     466             :             }
     467             :             else
     468             :             {
     469           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     470             :                          "R 'dim' dimension wrong.");
     471           0 :                 return nullptr;
     472             :             }
     473             :         }
     474           0 :         else if (nObjCode % 256 == R_REALSXP)
     475             :         {
     476           0 :             int nCount = poDS->ReadInteger();
     477           0 :             while (nCount > 0 && !VSIFEofL(poDS->fp))
     478             :             {
     479           0 :                 nCount--;
     480           0 :                 poDS->ReadFloat();
     481             :             }
     482             :         }
     483           0 :         else if (nObjCode % 256 == R_INTSXP)
     484             :         {
     485           0 :             int nCount = poDS->ReadInteger();
     486           0 :             while (nCount > 0 && !VSIFEofL(poDS->fp))
     487             :             {
     488           0 :                 nCount--;
     489           0 :                 poDS->ReadInteger();
     490             :             }
     491             :         }
     492           0 :         else if (nObjCode % 256 == R_STRSXP)
     493             :         {
     494           0 :             int nCount = poDS->ReadInteger();
     495           0 :             while (nCount > 0 && !VSIFEofL(poDS->fp))
     496             :             {
     497           0 :                 nCount--;
     498           0 :                 poDS->ReadString();
     499             :             }
     500             :         }
     501           0 :         else if (nObjCode % 256 == R_CHARSXP)
     502             :         {
     503           0 :             poDS->ReadString();
     504             :         }
     505             :     }
     506             : 
     507          10 :     if (poDS->nRasterXSize == 0)
     508             :     {
     509           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     510             :                  "Failed to find dim dimension information for R dataset.");
     511           0 :         return nullptr;
     512             :     }
     513             : 
     514          20 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
     515          10 :         !GDALCheckBandCount(nBandCount, TRUE))
     516             :     {
     517           0 :         return nullptr;
     518             :     }
     519             : 
     520          10 :     GIntBig result = 0;
     521          10 :     bool ok = SafeMult(nBandCount, poDS->nRasterXSize, &result);
     522          10 :     ok &= SafeMult(result, poDS->nRasterYSize, &result);
     523             : 
     524          10 :     if (!ok || nValueCount < result)
     525             :     {
     526           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Not enough pixel data.");
     527           0 :         return nullptr;
     528             :     }
     529             : 
     530             :     // Create the raster band object(s).
     531          24 :     for (int iBand = 0; iBand < nBandCount; iBand++)
     532             :     {
     533           0 :         std::unique_ptr<GDALRasterBand> poBand;
     534             : 
     535          14 :         if (poDS->bASCII)
     536           7 :             poBand = std::make_unique<RRasterBand>(
     537           7 :                 poDS.get(), iBand + 1,
     538          14 :                 poDS->padfMatrixValues + static_cast<size_t>(iBand) *
     539           7 :                                              poDS->nRasterXSize *
     540          14 :                                              poDS->nRasterYSize);
     541             :         else
     542             :         {
     543          14 :             poBand = RawRasterBand::Create(
     544           7 :                 poDS.get(), iBand + 1, poDS->fp,
     545          14 :                 poDS->nStartOfData +
     546           7 :                     static_cast<vsi_l_offset>(poDS->nRasterXSize) *
     547           7 :                         poDS->nRasterYSize * 8 * iBand,
     548           7 :                 8, poDS->nRasterXSize * 8, GDT_Float64,
     549             :                 RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN,
     550           7 :                 RawRasterBand::OwnFP::NO);
     551           7 :             if (!poBand)
     552           0 :                 return nullptr;
     553             :         }
     554          14 :         poDS->SetBand(iBand + 1, std::move(poBand));
     555             :     }
     556             : 
     557             :     // Initialize any PAM information.
     558          10 :     poDS->SetDescription(poOpenInfo->pszFilename);
     559          10 :     poDS->TryLoadXML();
     560             : 
     561             :     // Check for overviews.
     562          10 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     563             : 
     564          10 :     return poDS.release();
     565             : }
     566             : 
     567             : /************************************************************************/
     568             : /*                        GDALRegister_R()                              */
     569             : /************************************************************************/
     570             : 
     571        1686 : void GDALRegister_R()
     572             : 
     573             : {
     574        1686 :     if (GDALGetDriverByName("R") != nullptr)
     575         302 :         return;
     576             : 
     577        1384 :     GDALDriver *poDriver = new GDALDriver();
     578             : 
     579        1384 :     poDriver->SetDescription("R");
     580        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     581        1384 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "R Object Data Store");
     582        1384 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/r.html");
     583        1384 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "rda");
     584        1384 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");
     585        1384 :     poDriver->SetMetadataItem(
     586             :         GDAL_DMD_CREATIONOPTIONLIST,
     587             :         "<CreationOptionList>"
     588             :         "   <Option name='ASCII' type='boolean' description='For ASCII output, "
     589             :         "default NO'/>"
     590             :         "   <Option name='COMPRESS' type='boolean' description='Produced "
     591             :         "Compressed output, default YES'/>"
     592        1384 :         "</CreationOptionList>");
     593             : 
     594        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     595             : 
     596        1384 :     poDriver->pfnOpen = RDataset::Open;
     597        1384 :     poDriver->pfnIdentify = RDataset::Identify;
     598        1384 :     poDriver->pfnCreateCopy = RCreateCopy;
     599             : 
     600        1384 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     601             : }