1669 строки
58 KiB
C++
1669 строки
58 KiB
C++
/******************************************************************************
|
|
*
|
|
* Project: GDAL
|
|
* Purpose: Implements Arc/Info ASCII Grid Format.
|
|
* Author: Frank Warmerdam, warmerdam@pobox.com
|
|
*
|
|
******************************************************************************
|
|
* Copyright (c) 2001, Frank Warmerdam (warmerdam@pobox.com)
|
|
* Copyright (c) 2007-2012, Even Rouault <even dot rouault at spatialys.com>
|
|
* Copyright (c) 2014, Kyle Shannon <kyle at pobox dot com>
|
|
*
|
|
* Permission is hereby granted, free of charge, to any person obtaining a
|
|
* copy of this software and associated documentation files (the "Software"),
|
|
* to deal in the Software without restriction, including without limitation
|
|
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
|
* and/or sell copies of the Software, and to permit persons to whom the
|
|
* Software is furnished to do so, subject to the following conditions:
|
|
*
|
|
* The above copyright notice and this permission notice shall be included
|
|
* in all copies or substantial portions of the Software.
|
|
*
|
|
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
|
|
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
|
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
|
* DEALINGS IN THE SOFTWARE.
|
|
****************************************************************************/
|
|
|
|
// We need cpl_port as first include to avoid VSIStatBufL being not
|
|
// defined on i586-mingw32msvc.
|
|
#include "cpl_port.h"
|
|
#include "aaigriddataset.h"
|
|
#include "gdal_frmts.h"
|
|
|
|
#include <cassert>
|
|
#include <cctype>
|
|
#include <climits>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <cstdio>
|
|
#include <cstdlib>
|
|
#include <cstring>
|
|
#if HAVE_FCNTL_H
|
|
#include <fcntl.h>
|
|
#endif
|
|
|
|
#include <algorithm>
|
|
#include <limits>
|
|
#include <string>
|
|
|
|
#include "cpl_conv.h"
|
|
#include "cpl_error.h"
|
|
#include "cpl_progress.h"
|
|
#include "cpl_string.h"
|
|
#include "cpl_vsi.h"
|
|
#include "gdal.h"
|
|
#include "gdal_pam.h"
|
|
#include "gdal_priv.h"
|
|
#include "ogr_core.h"
|
|
#include "ogr_spatialref.h"
|
|
|
|
namespace
|
|
{
|
|
|
|
float DoubleToFloatClamp(double dfValue)
|
|
{
|
|
if (dfValue <= std::numeric_limits<float>::lowest())
|
|
return std::numeric_limits<float>::lowest();
|
|
if (dfValue >= std::numeric_limits<float>::max())
|
|
return std::numeric_limits<float>::max();
|
|
return static_cast<float>(dfValue);
|
|
}
|
|
|
|
// Cast to float and back for make sure the NoData value matches
|
|
// that expressed by a float value. Clamps to the range of a float
|
|
// if the value is too large. Preserves +/-inf and NaN.
|
|
// TODO(schwehr): This should probably be moved to port as it is likely
|
|
// to be needed for other formats.
|
|
double MapNoDataToFloat(double dfNoDataValue)
|
|
{
|
|
if (CPLIsInf(dfNoDataValue) || CPLIsNan(dfNoDataValue))
|
|
return dfNoDataValue;
|
|
|
|
if (dfNoDataValue >= std::numeric_limits<float>::max())
|
|
return std::numeric_limits<float>::max();
|
|
|
|
if (dfNoDataValue <= -std::numeric_limits<float>::max())
|
|
return -std::numeric_limits<float>::max();
|
|
|
|
return static_cast<double>(static_cast<float>(dfNoDataValue));
|
|
}
|
|
|
|
} // namespace
|
|
|
|
static CPLString OSR_GDS(char **papszNV, const char *pszField,
|
|
const char *pszDefaultValue);
|
|
|
|
/************************************************************************/
|
|
/* AAIGRasterBand() */
|
|
/************************************************************************/
|
|
|
|
AAIGRasterBand::AAIGRasterBand(AAIGDataset *poDSIn, int nDataStart)
|
|
: panLineOffset(nullptr)
|
|
{
|
|
poDS = poDSIn;
|
|
|
|
nBand = 1;
|
|
eDataType = poDSIn->eDataType;
|
|
|
|
nBlockXSize = poDSIn->nRasterXSize;
|
|
nBlockYSize = 1;
|
|
|
|
panLineOffset = static_cast<GUIntBig *>(
|
|
VSI_CALLOC_VERBOSE(poDSIn->nRasterYSize, sizeof(GUIntBig)));
|
|
if (panLineOffset == nullptr)
|
|
{
|
|
return;
|
|
}
|
|
panLineOffset[0] = nDataStart;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* ~AAIGRasterBand() */
|
|
/************************************************************************/
|
|
|
|
AAIGRasterBand::~AAIGRasterBand()
|
|
{
|
|
CPLFree(panLineOffset);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* IReadBlock() */
|
|
/************************************************************************/
|
|
|
|
CPLErr AAIGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
|
|
|
|
{
|
|
AAIGDataset *poODS = static_cast<AAIGDataset *>(poDS);
|
|
|
|
if (nBlockYOff < 0 || nBlockYOff > poODS->nRasterYSize - 1 ||
|
|
nBlockXOff != 0 || panLineOffset == nullptr || poODS->fp == nullptr)
|
|
return CE_Failure;
|
|
|
|
if (panLineOffset[nBlockYOff] == 0)
|
|
{
|
|
for (int iPrevLine = 1; iPrevLine <= nBlockYOff; iPrevLine++)
|
|
if (panLineOffset[iPrevLine] == 0)
|
|
IReadBlock(nBlockXOff, iPrevLine - 1, nullptr);
|
|
}
|
|
|
|
if (panLineOffset[nBlockYOff] == 0)
|
|
return CE_Failure;
|
|
|
|
if (poODS->Seek(panLineOffset[nBlockYOff]) != 0)
|
|
{
|
|
ReportError(CE_Failure, CPLE_FileIO,
|
|
"Can't seek to offset %lu in input file to read data.",
|
|
static_cast<long unsigned int>(panLineOffset[nBlockYOff]));
|
|
return CE_Failure;
|
|
}
|
|
|
|
for (int iPixel = 0; iPixel < poODS->nRasterXSize;)
|
|
{
|
|
// Suck up any pre-white space.
|
|
char chNext = '\0';
|
|
do
|
|
{
|
|
chNext = poODS->Getc();
|
|
} while (isspace(static_cast<unsigned char>(chNext)));
|
|
|
|
char szToken[500] = {'\0'};
|
|
int iTokenChar = 0;
|
|
while (chNext != '\0' && !isspace((unsigned char)chNext))
|
|
{
|
|
if (iTokenChar == sizeof(szToken) - 2)
|
|
{
|
|
ReportError(CE_Failure, CPLE_FileIO,
|
|
"Token too long at scanline %d.", nBlockYOff);
|
|
return CE_Failure;
|
|
}
|
|
|
|
szToken[iTokenChar++] = chNext;
|
|
chNext = poODS->Getc();
|
|
}
|
|
|
|
if (chNext == '\0' && (iPixel != poODS->nRasterXSize - 1 ||
|
|
nBlockYOff != poODS->nRasterYSize - 1))
|
|
{
|
|
ReportError(CE_Failure, CPLE_FileIO,
|
|
"File short, can't read line %d.", nBlockYOff);
|
|
return CE_Failure;
|
|
}
|
|
|
|
szToken[iTokenChar] = '\0';
|
|
|
|
if (pImage != nullptr)
|
|
{
|
|
// "null" seems to be specific of D12 software
|
|
// See https://github.com/OSGeo/gdal/issues/5095
|
|
if (eDataType == GDT_Float64)
|
|
{
|
|
if (strcmp(szToken, "null") == 0)
|
|
reinterpret_cast<double *>(pImage)[iPixel] =
|
|
-std::numeric_limits<double>::max();
|
|
else
|
|
reinterpret_cast<double *>(pImage)[iPixel] =
|
|
CPLAtofM(szToken);
|
|
}
|
|
else if (eDataType == GDT_Float32)
|
|
{
|
|
if (strcmp(szToken, "null") == 0)
|
|
reinterpret_cast<float *>(pImage)[iPixel] =
|
|
-std::numeric_limits<float>::max();
|
|
else
|
|
reinterpret_cast<float *>(pImage)[iPixel] =
|
|
DoubleToFloatClamp(CPLAtofM(szToken));
|
|
}
|
|
else
|
|
reinterpret_cast<GInt32 *>(pImage)[iPixel] =
|
|
static_cast<GInt32>(atoi(szToken));
|
|
}
|
|
|
|
iPixel++;
|
|
}
|
|
|
|
if (nBlockYOff < poODS->nRasterYSize - 1)
|
|
panLineOffset[nBlockYOff + 1] = poODS->Tell();
|
|
|
|
return CE_None;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GetNoDataValue() */
|
|
/************************************************************************/
|
|
|
|
double AAIGRasterBand::GetNoDataValue(int *pbSuccess)
|
|
|
|
{
|
|
AAIGDataset *poODS = static_cast<AAIGDataset *>(poDS);
|
|
|
|
if (pbSuccess)
|
|
*pbSuccess = poODS->bNoDataSet;
|
|
|
|
return poODS->dfNoDataValue;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* SetNoDataValue() */
|
|
/************************************************************************/
|
|
|
|
CPLErr AAIGRasterBand::SetNoDataValue(double dfNoData)
|
|
|
|
{
|
|
AAIGDataset *poODS = static_cast<AAIGDataset *>(poDS);
|
|
|
|
poODS->bNoDataSet = true;
|
|
poODS->dfNoDataValue = dfNoData;
|
|
|
|
return CE_None;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* ==================================================================== */
|
|
/* AAIGDataset */
|
|
/* ==================================================================== */
|
|
/************************************************************************/
|
|
|
|
/************************************************************************/
|
|
/* AAIGDataset() */
|
|
/************************************************************************/
|
|
|
|
AAIGDataset::AAIGDataset()
|
|
: fp(nullptr), papszPrj(nullptr), nBufferOffset(0), nOffsetInBuffer(256),
|
|
eDataType(GDT_Int32), bNoDataSet(false), dfNoDataValue(-9999.0)
|
|
{
|
|
m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
|
|
adfGeoTransform[0] = 0.0;
|
|
adfGeoTransform[1] = 1.0;
|
|
adfGeoTransform[2] = 0.0;
|
|
adfGeoTransform[3] = 0.0;
|
|
adfGeoTransform[4] = 0.0;
|
|
adfGeoTransform[5] = 1.0;
|
|
memset(achReadBuf, 0, sizeof(achReadBuf));
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* ~AAIGDataset() */
|
|
/************************************************************************/
|
|
|
|
AAIGDataset::~AAIGDataset()
|
|
|
|
{
|
|
FlushCache(true);
|
|
|
|
if (fp != nullptr)
|
|
{
|
|
if (VSIFCloseL(fp) != 0)
|
|
{
|
|
ReportError(CE_Failure, CPLE_FileIO, "I/O error");
|
|
}
|
|
}
|
|
|
|
CSLDestroy(papszPrj);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Tell() */
|
|
/************************************************************************/
|
|
|
|
GUIntBig AAIGDataset::Tell() const
|
|
{
|
|
return nBufferOffset + nOffsetInBuffer;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Seek() */
|
|
/************************************************************************/
|
|
|
|
int AAIGDataset::Seek(GUIntBig nNewOffset)
|
|
|
|
{
|
|
nOffsetInBuffer = sizeof(achReadBuf);
|
|
return VSIFSeekL(fp, nNewOffset, SEEK_SET);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Getc() */
|
|
/* */
|
|
/* Read a single character from the input file (efficiently we */
|
|
/* hope). */
|
|
/************************************************************************/
|
|
|
|
char AAIGDataset::Getc()
|
|
|
|
{
|
|
if (nOffsetInBuffer < static_cast<int>(sizeof(achReadBuf)))
|
|
return achReadBuf[nOffsetInBuffer++];
|
|
|
|
nBufferOffset = VSIFTellL(fp);
|
|
const int nRead =
|
|
static_cast<int>(VSIFReadL(achReadBuf, 1, sizeof(achReadBuf), fp));
|
|
for (unsigned int i = nRead; i < sizeof(achReadBuf); i++)
|
|
achReadBuf[i] = '\0';
|
|
|
|
nOffsetInBuffer = 0;
|
|
|
|
return achReadBuf[nOffsetInBuffer++];
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GetFileList() */
|
|
/************************************************************************/
|
|
|
|
char **AAIGDataset::GetFileList()
|
|
|
|
{
|
|
char **papszFileList = GDALPamDataset::GetFileList();
|
|
|
|
if (papszPrj != nullptr)
|
|
papszFileList = CSLAddString(papszFileList, osPrjFilename);
|
|
|
|
return papszFileList;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Identify() */
|
|
/************************************************************************/
|
|
|
|
int AAIGDataset::Identify(GDALOpenInfo *poOpenInfo)
|
|
|
|
{
|
|
// Does this look like an AI grid file?
|
|
if (poOpenInfo->nHeaderBytes < 40 ||
|
|
!(STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "ncols") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "nrows") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "xllcorner") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "yllcorner") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "xllcenter") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "yllcenter") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "dx") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "dy") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "cellsize")))
|
|
return FALSE;
|
|
|
|
return TRUE;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Identify() */
|
|
/************************************************************************/
|
|
|
|
int GRASSASCIIDataset::Identify(GDALOpenInfo *poOpenInfo)
|
|
|
|
{
|
|
// Does this look like a GRASS ASCII grid file?
|
|
if (poOpenInfo->nHeaderBytes < 40 ||
|
|
!(STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "north:") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "south:") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "east:") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "west:") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "rows:") ||
|
|
STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "cols:")))
|
|
return FALSE;
|
|
|
|
return TRUE;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Identify() */
|
|
/************************************************************************/
|
|
|
|
int ISGDataset::Identify(GDALOpenInfo *poOpenInfo)
|
|
|
|
{
|
|
// Does this look like a ISG grid file?
|
|
if (poOpenInfo->nHeaderBytes < 40 ||
|
|
!(strstr((const char *)poOpenInfo->pabyHeader, "model name") !=
|
|
nullptr &&
|
|
strstr((const char *)poOpenInfo->pabyHeader, "lat min") != nullptr &&
|
|
strstr((const char *)poOpenInfo->pabyHeader, "lat max") != nullptr &&
|
|
strstr((const char *)poOpenInfo->pabyHeader, "lon min") != nullptr &&
|
|
strstr((const char *)poOpenInfo->pabyHeader, "lon max") != nullptr &&
|
|
strstr((const char *)poOpenInfo->pabyHeader, "nrows") != nullptr &&
|
|
strstr((const char *)poOpenInfo->pabyHeader, "ncols") != nullptr))
|
|
return FALSE;
|
|
|
|
return TRUE;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Open() */
|
|
/************************************************************************/
|
|
|
|
GDALDataset *AAIGDataset::Open(GDALOpenInfo *poOpenInfo)
|
|
{
|
|
#ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
|
|
// During fuzzing, do not use Identify to reject crazy content.
|
|
if (!Identify(poOpenInfo))
|
|
return nullptr;
|
|
#endif
|
|
|
|
return CommonOpen(poOpenInfo, FORMAT_AAIG);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* ParseHeader() */
|
|
/************************************************************************/
|
|
|
|
int AAIGDataset::ParseHeader(const char *pszHeader, const char *pszDataType)
|
|
{
|
|
char **papszTokens = CSLTokenizeString2(pszHeader, " \n\r\t", 0);
|
|
const int nTokens = CSLCount(papszTokens);
|
|
|
|
int i = 0;
|
|
if ((i = CSLFindString(papszTokens, "ncols")) < 0 || i + 1 >= nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
nRasterXSize = atoi(papszTokens[i + 1]);
|
|
if ((i = CSLFindString(papszTokens, "nrows")) < 0 || i + 1 >= nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
nRasterYSize = atoi(papszTokens[i + 1]);
|
|
|
|
if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
|
|
// TODO(schwehr): Would be good to also factor the file size into the max.
|
|
// TODO(schwehr): Allow the user to disable this check.
|
|
// The driver allocates a panLineOffset array based on nRasterYSize
|
|
constexpr int kMaxDimSize = 10000000; // 1e7 cells.
|
|
if (nRasterXSize > kMaxDimSize || nRasterYSize > kMaxDimSize)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
|
|
double dfCellDX = 0.0;
|
|
double dfCellDY = 0.0;
|
|
if ((i = CSLFindString(papszTokens, "cellsize")) < 0)
|
|
{
|
|
int iDX, iDY;
|
|
if ((iDX = CSLFindString(papszTokens, "dx")) < 0 ||
|
|
(iDY = CSLFindString(papszTokens, "dy")) < 0 ||
|
|
iDX + 1 >= nTokens || iDY + 1 >= nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
|
|
dfCellDX = CPLAtofM(papszTokens[iDX + 1]);
|
|
dfCellDY = CPLAtofM(papszTokens[iDY + 1]);
|
|
}
|
|
else
|
|
{
|
|
if (i + 1 >= nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
dfCellDY = CPLAtofM(papszTokens[i + 1]);
|
|
dfCellDX = dfCellDY;
|
|
}
|
|
|
|
int j = 0;
|
|
if ((i = CSLFindString(papszTokens, "xllcorner")) >= 0 &&
|
|
(j = CSLFindString(papszTokens, "yllcorner")) >= 0 && i + 1 < nTokens &&
|
|
j + 1 < nTokens)
|
|
{
|
|
adfGeoTransform[0] = CPLAtofM(papszTokens[i + 1]);
|
|
|
|
// Small hack to compensate from insufficient precision in cellsize
|
|
// parameter in datasets of
|
|
// http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3
|
|
if ((nRasterXSize % 360) == 0 &&
|
|
fabs(adfGeoTransform[0] - (-180.0)) < 1e-12 &&
|
|
dfCellDX == dfCellDY &&
|
|
fabs(dfCellDX - (360.0 / nRasterXSize)) < 1e-9)
|
|
{
|
|
dfCellDY = 360.0 / nRasterXSize;
|
|
dfCellDX = dfCellDY;
|
|
}
|
|
|
|
adfGeoTransform[1] = dfCellDX;
|
|
adfGeoTransform[2] = 0.0;
|
|
adfGeoTransform[3] =
|
|
CPLAtofM(papszTokens[j + 1]) + nRasterYSize * dfCellDY;
|
|
adfGeoTransform[4] = 0.0;
|
|
adfGeoTransform[5] = -dfCellDY;
|
|
}
|
|
else if ((i = CSLFindString(papszTokens, "xllcenter")) >= 0 &&
|
|
(j = CSLFindString(papszTokens, "yllcenter")) >= 0 &&
|
|
i + 1 < nTokens && j + 1 < nTokens)
|
|
{
|
|
SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
|
|
|
|
adfGeoTransform[0] = CPLAtofM(papszTokens[i + 1]) - 0.5 * dfCellDX;
|
|
adfGeoTransform[1] = dfCellDX;
|
|
adfGeoTransform[2] = 0.0;
|
|
adfGeoTransform[3] = CPLAtofM(papszTokens[j + 1]) - 0.5 * dfCellDY +
|
|
nRasterYSize * dfCellDY;
|
|
adfGeoTransform[4] = 0.0;
|
|
adfGeoTransform[5] = -dfCellDY;
|
|
}
|
|
else
|
|
{
|
|
adfGeoTransform[0] = 0.0;
|
|
adfGeoTransform[1] = dfCellDX;
|
|
adfGeoTransform[2] = 0.0;
|
|
adfGeoTransform[3] = 0.0;
|
|
adfGeoTransform[4] = 0.0;
|
|
adfGeoTransform[5] = -dfCellDY;
|
|
}
|
|
|
|
if ((i = CSLFindString(papszTokens, "NODATA_value")) >= 0 &&
|
|
i + 1 < nTokens)
|
|
{
|
|
const char *pszNoData = papszTokens[i + 1];
|
|
|
|
bNoDataSet = true;
|
|
if (strcmp(pszNoData, "null") == 0)
|
|
{
|
|
// "null" seems to be specific of D12 software
|
|
// See https://github.com/OSGeo/gdal/issues/5095
|
|
if (pszDataType == nullptr || eDataType == GDT_Float32)
|
|
{
|
|
dfNoDataValue = -std::numeric_limits<float>::max();
|
|
eDataType = GDT_Float32;
|
|
}
|
|
else
|
|
{
|
|
dfNoDataValue = -std::numeric_limits<double>::max();
|
|
eDataType = GDT_Float64;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
dfNoDataValue = CPLAtofM(pszNoData);
|
|
if (pszDataType == nullptr &&
|
|
(strchr(pszNoData, '.') != nullptr ||
|
|
strchr(pszNoData, ',') != nullptr ||
|
|
std::numeric_limits<int>::min() > dfNoDataValue ||
|
|
dfNoDataValue > std::numeric_limits<int>::max()))
|
|
{
|
|
eDataType = GDT_Float32;
|
|
if (!CPLIsInf(dfNoDataValue) &&
|
|
(fabs(dfNoDataValue) < std::numeric_limits<float>::min() ||
|
|
fabs(dfNoDataValue) > std::numeric_limits<float>::max()))
|
|
{
|
|
eDataType = GDT_Float64;
|
|
}
|
|
}
|
|
if (eDataType == GDT_Float32)
|
|
{
|
|
dfNoDataValue = MapNoDataToFloat(dfNoDataValue);
|
|
}
|
|
}
|
|
}
|
|
|
|
CSLDestroy(papszTokens);
|
|
|
|
return TRUE;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Open() */
|
|
/************************************************************************/
|
|
|
|
GDALDataset *GRASSASCIIDataset::Open(GDALOpenInfo *poOpenInfo)
|
|
{
|
|
#ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
|
|
// During fuzzing, do not use Identify to reject crazy content.
|
|
if (!Identify(poOpenInfo))
|
|
return nullptr;
|
|
#endif
|
|
|
|
return CommonOpen(poOpenInfo, FORMAT_GRASSASCII);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* ParseHeader() */
|
|
/************************************************************************/
|
|
|
|
int GRASSASCIIDataset::ParseHeader(const char *pszHeader,
|
|
const char *pszDataType)
|
|
{
|
|
char **papszTokens = CSLTokenizeString2(pszHeader, " \n\r\t:", 0);
|
|
const int nTokens = CSLCount(papszTokens);
|
|
int i = 0;
|
|
if ((i = CSLFindString(papszTokens, "cols")) < 0 || i + 1 >= nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
nRasterXSize = atoi(papszTokens[i + 1]);
|
|
if ((i = CSLFindString(papszTokens, "rows")) < 0 || i + 1 >= nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
nRasterYSize = atoi(papszTokens[i + 1]);
|
|
|
|
if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
|
|
// TODO(schwehr): Would be good to also factor the file size into the max.
|
|
// TODO(schwehr): Allow the user to disable this check.
|
|
// The driver allocates a panLineOffset array based on nRasterYSize
|
|
constexpr int kMaxDimSize = 10000000; // 1e7 cells.
|
|
if (nRasterXSize > kMaxDimSize || nRasterYSize > kMaxDimSize)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
|
|
const int iNorth = CSLFindString(papszTokens, "north");
|
|
const int iSouth = CSLFindString(papszTokens, "south");
|
|
const int iEast = CSLFindString(papszTokens, "east");
|
|
const int iWest = CSLFindString(papszTokens, "west");
|
|
|
|
if (iNorth == -1 || iSouth == -1 || iEast == -1 || iWest == -1 ||
|
|
std::max(std::max(iNorth, iSouth), std::max(iEast, iWest)) + 1 >=
|
|
nTokens)
|
|
{
|
|
CSLDestroy(papszTokens);
|
|
return FALSE;
|
|
}
|
|
|
|
const double dfNorth = CPLAtofM(papszTokens[iNorth + 1]);
|
|
const double dfSouth = CPLAtofM(papszTokens[iSouth + 1]);
|
|
const double dfEast = CPLAtofM(papszTokens[iEast + 1]);
|
|
const double dfWest = CPLAtofM(papszTokens[iWest + 1]);
|
|
const double dfPixelXSize = (dfEast - dfWest) / nRasterXSize;
|
|
const double dfPixelYSize = (dfNorth - dfSouth) / nRasterYSize;
|
|
|
|
adfGeoTransform[0] = dfWest;
|
|
adfGeoTransform[1] = dfPixelXSize;
|
|
adfGeoTransform[2] = 0.0;
|
|
adfGeoTransform[3] = dfNorth;
|
|
adfGeoTransform[4] = 0.0;
|
|
adfGeoTransform[5] = -dfPixelYSize;
|
|
|
|
if ((i = CSLFindString(papszTokens, "null")) >= 0 && i + 1 < nTokens)
|
|
{
|
|
const char *pszNoData = papszTokens[i + 1];
|
|
|
|
bNoDataSet = true;
|
|
dfNoDataValue = CPLAtofM(pszNoData);
|
|
if (pszDataType == nullptr &&
|
|
(strchr(pszNoData, '.') != nullptr ||
|
|
strchr(pszNoData, ',') != nullptr ||
|
|
std::numeric_limits<int>::min() > dfNoDataValue ||
|
|
dfNoDataValue > std::numeric_limits<int>::max()))
|
|
{
|
|
eDataType = GDT_Float32;
|
|
}
|
|
if (eDataType == GDT_Float32)
|
|
{
|
|
dfNoDataValue = MapNoDataToFloat(dfNoDataValue);
|
|
}
|
|
}
|
|
|
|
if ((i = CSLFindString(papszTokens, "type")) >= 0 && i + 1 < nTokens)
|
|
{
|
|
const char *pszType = papszTokens[i + 1];
|
|
if (EQUAL(pszType, "int"))
|
|
eDataType = GDT_Int32;
|
|
else if (EQUAL(pszType, "float"))
|
|
eDataType = GDT_Float32;
|
|
else if (EQUAL(pszType, "double"))
|
|
eDataType = GDT_Float64;
|
|
else
|
|
{
|
|
ReportError(CE_Warning, CPLE_AppDefined,
|
|
"Invalid value for type parameter : %s", pszType);
|
|
}
|
|
}
|
|
|
|
CSLDestroy(papszTokens);
|
|
|
|
return TRUE;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* Open() */
|
|
/************************************************************************/
|
|
|
|
GDALDataset *ISGDataset::Open(GDALOpenInfo *poOpenInfo)
|
|
{
|
|
#ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
|
|
// During fuzzing, do not use Identify to reject crazy content.
|
|
if (!Identify(poOpenInfo))
|
|
return nullptr;
|
|
#endif
|
|
|
|
return CommonOpen(poOpenInfo, FORMAT_ISG);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* ParseHeader() */
|
|
/************************************************************************/
|
|
|
|
int ISGDataset::ParseHeader(const char *pszHeader, const char *)
|
|
{
|
|
// See http://www.isgeoid.polimi.it/Geoid/ISG_format_20160121.pdf
|
|
|
|
CPLStringList aosLines(CSLTokenizeString2(pszHeader, "\n\r", 0));
|
|
CPLString osLatMin;
|
|
CPLString osLatMax;
|
|
CPLString osLonMin;
|
|
CPLString osLonMax;
|
|
CPLString osDeltaLat;
|
|
CPLString osDeltaLon;
|
|
CPLString osRows;
|
|
CPLString osCols;
|
|
CPLString osNodata;
|
|
for (int iLine = 0; iLine < aosLines.size(); iLine++)
|
|
{
|
|
CPLStringList aosTokens(CSLTokenizeString2(aosLines[iLine], ":=", 0));
|
|
if (aosTokens.size() == 2)
|
|
{
|
|
CPLString osLeft(aosTokens[0]);
|
|
osLeft.Trim();
|
|
CPLString osRight(aosTokens[1]);
|
|
osRight.Trim();
|
|
if (osLeft == "lat min")
|
|
osLatMin = osRight;
|
|
else if (osLeft == "lat max")
|
|
osLatMax = osRight;
|
|
else if (osLeft == "lon min")
|
|
osLonMin = osRight;
|
|
else if (osLeft == "lon max")
|
|
osLonMax = osRight;
|
|
else if (osLeft == "delta lat")
|
|
osDeltaLat = osRight;
|
|
else if (osLeft == "delta lon")
|
|
osDeltaLon = osRight;
|
|
else if (osLeft == "nrows")
|
|
osRows = osRight;
|
|
else if (osLeft == "ncols")
|
|
osCols = osRight;
|
|
else if (osLeft == "nodata")
|
|
osNodata = osRight;
|
|
else if (osLeft == "model name")
|
|
SetMetadataItem("MODEL_NAME", osRight);
|
|
else if (osLeft == "model type")
|
|
SetMetadataItem("MODEL_TYPE", osRight);
|
|
else if (osLeft == "units")
|
|
osUnits = osRight;
|
|
}
|
|
}
|
|
if (osLatMin.empty() || osLatMax.empty() || osLonMin.empty() ||
|
|
osLonMax.empty() || osDeltaLat.empty() || osDeltaLon.empty() ||
|
|
osRows.empty() || osCols.empty())
|
|
{
|
|
return FALSE;
|
|
}
|
|
double dfLatMin = CPLAtof(osLatMin);
|
|
double dfLatMax = CPLAtof(osLatMax);
|
|
double dfLonMin = CPLAtof(osLonMin);
|
|
double dfLonMax = CPLAtof(osLonMax);
|
|
double dfDeltaLon = CPLAtof(osDeltaLon);
|
|
double dfDeltaLat = CPLAtof(osDeltaLat);
|
|
const int nRows = atoi(osRows);
|
|
const int nCols = atoi(osCols);
|
|
if (nRows <= 0 || nCols <= 0 ||
|
|
!(dfDeltaLat > 0 && dfDeltaLon > 0 && dfDeltaLat < 180 &&
|
|
dfDeltaLon < 360))
|
|
{
|
|
return FALSE;
|
|
}
|
|
|
|
// Correct rounding errors.
|
|
|
|
const auto TryRoundTo = [](double &dfDelta, double dfRoundedDelta,
|
|
double &dfMin, double &dfMax, int nVals,
|
|
double dfRelTol)
|
|
{
|
|
double dfMinTry = dfMin;
|
|
double dfMaxTry = dfMax;
|
|
double dfDeltaTry = dfDelta;
|
|
if (dfRoundedDelta != dfDelta &&
|
|
fabs(fabs(dfMin / dfRoundedDelta) -
|
|
(floor(fabs(dfMin / dfRoundedDelta)) + 0.5)) < dfRelTol &&
|
|
fabs(fabs(dfMax / dfRoundedDelta) -
|
|
(floor(fabs(dfMax / dfRoundedDelta)) + 0.5)) < dfRelTol)
|
|
{
|
|
{
|
|
double dfVal = (floor(fabs(dfMin / dfRoundedDelta)) + 0.5) *
|
|
dfRoundedDelta;
|
|
dfMinTry = (dfMin < 0) ? -dfVal : dfVal;
|
|
}
|
|
{
|
|
double dfVal = (floor(fabs(dfMax / dfRoundedDelta)) + 0.5) *
|
|
dfRoundedDelta;
|
|
dfMaxTry = (dfMax < 0) ? -dfVal : dfVal;
|
|
}
|
|
dfDeltaTry = dfRoundedDelta;
|
|
}
|
|
else if (dfRoundedDelta != dfDelta &&
|
|
fabs(fabs(dfMin / dfRoundedDelta) -
|
|
(floor(fabs(dfMin / dfRoundedDelta) + 0.5) + 0.)) <
|
|
dfRelTol &&
|
|
fabs(fabs(dfMax / dfRoundedDelta) -
|
|
(floor(fabs(dfMax / dfRoundedDelta) + 0.5) + 0.)) <
|
|
dfRelTol)
|
|
{
|
|
{
|
|
double dfVal =
|
|
(floor(fabs(dfMin / dfRoundedDelta) + 0.5) + 0.) *
|
|
dfRoundedDelta;
|
|
dfMinTry = (dfMin < 0) ? -dfVal : dfVal;
|
|
}
|
|
{
|
|
double dfVal =
|
|
(floor(fabs(dfMax / dfRoundedDelta) + 0.5) + 0.) *
|
|
dfRoundedDelta;
|
|
dfMaxTry = (dfMax < 0) ? -dfVal : dfVal;
|
|
}
|
|
dfDeltaTry = dfRoundedDelta;
|
|
}
|
|
if (fabs(dfMinTry + dfDeltaTry * nVals - dfMaxTry) <
|
|
dfRelTol * dfDeltaTry)
|
|
{
|
|
dfMin = dfMinTry;
|
|
dfMax = dfMaxTry;
|
|
dfDelta = dfDeltaTry;
|
|
return true;
|
|
}
|
|
return false;
|
|
};
|
|
|
|
const double dfRoundedDeltaLon =
|
|
(osDeltaLon == "0.0167" ||
|
|
(dfDeltaLon < 1 &&
|
|
fabs(1. / dfDeltaLon - floor(1. / dfDeltaLon + 0.5)) < 0.06))
|
|
? 1. / floor(1. / dfDeltaLon + 0.5)
|
|
: dfDeltaLon;
|
|
|
|
const double dfRoundedDeltaLat =
|
|
(osDeltaLat == "0.0167" ||
|
|
(dfDeltaLat < 1 &&
|
|
fabs(1. / dfDeltaLat - floor(1. / dfDeltaLat + 0.5)) < 0.06))
|
|
? 1. / floor(1. / dfDeltaLat + 0.5)
|
|
: dfDeltaLat;
|
|
|
|
bool bOK = TryRoundTo(dfDeltaLon, dfRoundedDeltaLon, dfLonMin, dfLonMax,
|
|
nCols, 1e-2) &&
|
|
TryRoundTo(dfDeltaLat, dfRoundedDeltaLat, dfLatMin, dfLatMax,
|
|
nRows, 1e-2);
|
|
if (!bOK && osDeltaLon == "0.0167" && osDeltaLat == "0.0167")
|
|
{
|
|
// For https://www.isgeoid.polimi.it/Geoid/America/Argentina/public/GEOIDEAR16_20160419.isg
|
|
bOK =
|
|
TryRoundTo(dfDeltaLon, 0.016667, dfLonMin, dfLonMax, nCols, 1e-1) &&
|
|
TryRoundTo(dfDeltaLat, 0.016667, dfLatMin, dfLatMax, nRows, 1e-1);
|
|
}
|
|
if (!bOK)
|
|
{
|
|
// 0.005 is what would be needed for the above GEOIDEAR16_20160419.isg
|
|
// file without the specific fine tuning done.
|
|
if ((fabs((dfLonMax - dfLonMin) / nCols - dfDeltaLon) <
|
|
0.005 * dfDeltaLon &&
|
|
fabs((dfLatMax - dfLatMin) / nRows - dfDeltaLat) <
|
|
0.005 * dfDeltaLat) ||
|
|
CPLTestBool(
|
|
CPLGetConfigOption("ISG_SKIP_GEOREF_CONSISTENCY_CHECK", "NO")))
|
|
{
|
|
CPLError(CE_Warning, CPLE_AppDefined,
|
|
"Georeference might be slightly approximate due to "
|
|
"rounding of coordinates and resolution in file header.");
|
|
dfDeltaLon = (dfLonMax - dfLonMin) / nCols;
|
|
dfDeltaLat = (dfLatMax - dfLatMin) / nRows;
|
|
}
|
|
else
|
|
{
|
|
CPLError(CE_Failure, CPLE_AppDefined,
|
|
"Inconsistent extent/resolution/raster dimension, or "
|
|
"rounding of coordinates and resolution in file header "
|
|
"higher than accepted. You may skip this consistency "
|
|
"check by setting the ISG_SKIP_GEOREF_CONSISTENCY_CHECK "
|
|
"configuration option to YES.");
|
|
return false;
|
|
}
|
|
}
|
|
nRasterXSize = nCols;
|
|
nRasterYSize = nRows;
|
|
adfGeoTransform[0] = dfLonMin;
|
|
adfGeoTransform[1] = dfDeltaLon;
|
|
adfGeoTransform[2] = 0.0;
|
|
adfGeoTransform[3] = dfLatMax;
|
|
adfGeoTransform[4] = 0.0;
|
|
adfGeoTransform[5] = -dfDeltaLat;
|
|
if (!osNodata.empty())
|
|
{
|
|
bNoDataSet = true;
|
|
dfNoDataValue = MapNoDataToFloat(CPLAtof(osNodata));
|
|
}
|
|
return TRUE;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* CommonOpen() */
|
|
/************************************************************************/
|
|
|
|
GDALDataset *AAIGDataset::CommonOpen(GDALOpenInfo *poOpenInfo,
|
|
GridFormat eFormat)
|
|
{
|
|
if (poOpenInfo->fpL == nullptr)
|
|
return nullptr;
|
|
|
|
// Create a corresponding GDALDataset.
|
|
AAIGDataset *poDS = nullptr;
|
|
|
|
if (eFormat == FORMAT_AAIG)
|
|
poDS = new AAIGDataset();
|
|
else if (eFormat == FORMAT_GRASSASCII)
|
|
poDS = new GRASSASCIIDataset();
|
|
else
|
|
{
|
|
poDS = new ISGDataset();
|
|
poDS->eDataType = GDT_Float32;
|
|
}
|
|
|
|
const char *pszDataTypeOption = eFormat == FORMAT_AAIG ? "AAIGRID_DATATYPE"
|
|
: eFormat == FORMAT_GRASSASCII
|
|
? "GRASSASCIIGRID_DATATYPE"
|
|
: nullptr;
|
|
|
|
const char *pszDataType =
|
|
pszDataTypeOption ? CPLGetConfigOption(pszDataTypeOption, nullptr)
|
|
: nullptr;
|
|
if (pszDataType == nullptr)
|
|
{
|
|
pszDataType =
|
|
CSLFetchNameValue(poOpenInfo->papszOpenOptions, "DATATYPE");
|
|
}
|
|
if (pszDataType != nullptr)
|
|
{
|
|
poDS->eDataType = GDALGetDataTypeByName(pszDataType);
|
|
if (!(poDS->eDataType == GDT_Int32 || poDS->eDataType == GDT_Float32 ||
|
|
poDS->eDataType == GDT_Float64))
|
|
{
|
|
ReportError(poOpenInfo->pszFilename, CE_Warning, CPLE_NotSupported,
|
|
"Unsupported value for %s : %s", pszDataTypeOption,
|
|
pszDataType);
|
|
poDS->eDataType = GDT_Int32;
|
|
pszDataType = nullptr;
|
|
}
|
|
}
|
|
|
|
// Parse the header.
|
|
if (!poDS->ParseHeader((const char *)poOpenInfo->pabyHeader, pszDataType))
|
|
{
|
|
delete poDS;
|
|
return nullptr;
|
|
}
|
|
|
|
poDS->fp = poOpenInfo->fpL;
|
|
poOpenInfo->fpL = nullptr;
|
|
|
|
// Find the start of real data.
|
|
int nStartOfData = 0;
|
|
|
|
if (eFormat == FORMAT_ISG)
|
|
{
|
|
const char *pszEOH =
|
|
strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
|
|
"end_of_head");
|
|
if (pszEOH == nullptr)
|
|
{
|
|
delete poDS;
|
|
return nullptr;
|
|
}
|
|
for (int i = 0; pszEOH[i]; i++)
|
|
{
|
|
if (pszEOH[i] == '\n' || pszEOH[i] == '\r')
|
|
{
|
|
nStartOfData =
|
|
static_cast<int>(pszEOH - reinterpret_cast<const char *>(
|
|
poOpenInfo->pabyHeader)) +
|
|
i;
|
|
break;
|
|
}
|
|
}
|
|
if (nStartOfData == 0)
|
|
{
|
|
delete poDS;
|
|
return nullptr;
|
|
}
|
|
if (poOpenInfo->pabyHeader[nStartOfData] == '\n' ||
|
|
poOpenInfo->pabyHeader[nStartOfData] == '\r')
|
|
{
|
|
nStartOfData++;
|
|
}
|
|
|
|
poDS->m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
|
|
}
|
|
else
|
|
{
|
|
for (int i = 2; true; i++)
|
|
{
|
|
if (poOpenInfo->pabyHeader[i] == '\0')
|
|
{
|
|
ReportError(poOpenInfo->pszFilename, CE_Failure,
|
|
CPLE_AppDefined,
|
|
"Couldn't find data values in ASCII Grid file.");
|
|
delete poDS;
|
|
return nullptr;
|
|
}
|
|
|
|
if (poOpenInfo->pabyHeader[i - 1] == '\n' ||
|
|
poOpenInfo->pabyHeader[i - 2] == '\n' ||
|
|
poOpenInfo->pabyHeader[i - 1] == '\r' ||
|
|
poOpenInfo->pabyHeader[i - 2] == '\r')
|
|
{
|
|
if ((!isalpha(poOpenInfo->pabyHeader[i]) ||
|
|
// null seems to be specific of D12 software
|
|
// See https://github.com/OSGeo/gdal/issues/5095
|
|
(i + 5 < poOpenInfo->nHeaderBytes &&
|
|
memcmp(poOpenInfo->pabyHeader + i, "null ", 5) == 0)) &&
|
|
poOpenInfo->pabyHeader[i] != '\n' &&
|
|
poOpenInfo->pabyHeader[i] != '\r')
|
|
{
|
|
nStartOfData = i;
|
|
|
|
// Beginning of real data found.
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Recognize the type of data.
|
|
CPLAssert(nullptr != poDS->fp);
|
|
|
|
if (pszDataType == nullptr && poDS->eDataType != GDT_Float32 &&
|
|
poDS->eDataType != GDT_Float64)
|
|
{
|
|
// Allocate 100K chunk + 1 extra byte for NULL character.
|
|
constexpr size_t nChunkSize = 1024 * 100;
|
|
GByte *pabyChunk = static_cast<GByte *>(
|
|
VSI_CALLOC_VERBOSE(nChunkSize + 1, sizeof(GByte)));
|
|
if (pabyChunk == nullptr)
|
|
{
|
|
delete poDS;
|
|
return nullptr;
|
|
}
|
|
pabyChunk[nChunkSize] = '\0';
|
|
|
|
if (VSIFSeekL(poDS->fp, nStartOfData, SEEK_SET) < 0)
|
|
{
|
|
delete poDS;
|
|
VSIFree(pabyChunk);
|
|
return nullptr;
|
|
}
|
|
|
|
// Scan for dot in subsequent chunks of data.
|
|
while (!VSIFEofL(poDS->fp))
|
|
{
|
|
const size_t nLen = VSIFReadL(pabyChunk, 1, nChunkSize, poDS->fp);
|
|
|
|
for (size_t i = 0; i < nLen; i++)
|
|
{
|
|
const GByte ch = pabyChunk[i];
|
|
if (ch == '.' || ch == ',' || ch == 'e' || ch == 'E')
|
|
{
|
|
poDS->eDataType = GDT_Float32;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Deallocate chunk.
|
|
VSIFree(pabyChunk);
|
|
}
|
|
|
|
// Create band information objects.
|
|
AAIGRasterBand *band = new AAIGRasterBand(poDS, nStartOfData);
|
|
poDS->SetBand(1, band);
|
|
if (band->panLineOffset == nullptr)
|
|
{
|
|
delete poDS;
|
|
return nullptr;
|
|
}
|
|
if (!poDS->osUnits.empty())
|
|
{
|
|
poDS->GetRasterBand(1)->SetUnitType(poDS->osUnits);
|
|
}
|
|
|
|
// Try to read projection file.
|
|
char *const pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
|
|
char *const pszBasename =
|
|
CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
|
|
|
|
poDS->osPrjFilename = CPLFormFilename(pszDirname, pszBasename, "prj");
|
|
int nRet = 0;
|
|
{
|
|
VSIStatBufL sStatBuf;
|
|
nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
|
|
}
|
|
if (nRet != 0 && VSIIsCaseSensitiveFS(poDS->osPrjFilename))
|
|
{
|
|
poDS->osPrjFilename = CPLFormFilename(pszDirname, pszBasename, "PRJ");
|
|
|
|
VSIStatBufL sStatBuf;
|
|
nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
|
|
}
|
|
|
|
if (nRet == 0)
|
|
{
|
|
poDS->papszPrj = CSLLoad(poDS->osPrjFilename);
|
|
|
|
CPLDebug("AAIGrid", "Loaded SRS from %s", poDS->osPrjFilename.c_str());
|
|
|
|
OGRSpatialReference oSRS;
|
|
oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
|
|
if (oSRS.importFromESRI(poDS->papszPrj) == OGRERR_NONE)
|
|
{
|
|
// If geographic values are in seconds, we must transform.
|
|
// Is there a code for minutes too?
|
|
if (oSRS.IsGeographic() &&
|
|
EQUAL(OSR_GDS(poDS->papszPrj, "Units", ""), "DS"))
|
|
{
|
|
poDS->adfGeoTransform[0] /= 3600.0;
|
|
poDS->adfGeoTransform[1] /= 3600.0;
|
|
poDS->adfGeoTransform[2] /= 3600.0;
|
|
poDS->adfGeoTransform[3] /= 3600.0;
|
|
poDS->adfGeoTransform[4] /= 3600.0;
|
|
poDS->adfGeoTransform[5] /= 3600.0;
|
|
}
|
|
|
|
poDS->m_oSRS = oSRS;
|
|
}
|
|
}
|
|
|
|
CPLFree(pszDirname);
|
|
CPLFree(pszBasename);
|
|
|
|
// Initialize any PAM information.
|
|
poDS->SetDescription(poOpenInfo->pszFilename);
|
|
poDS->TryLoadXML();
|
|
|
|
// Check for external overviews.
|
|
poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
|
|
poOpenInfo->GetSiblingFiles());
|
|
|
|
return poDS;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GetGeoTransform() */
|
|
/************************************************************************/
|
|
|
|
CPLErr AAIGDataset::GetGeoTransform(double *padfTransform)
|
|
|
|
{
|
|
memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
|
|
return CE_None;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GetSpatialRef() */
|
|
/************************************************************************/
|
|
|
|
const OGRSpatialReference *AAIGDataset::GetSpatialRef() const
|
|
{
|
|
return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* CreateCopy() */
|
|
/************************************************************************/
|
|
|
|
GDALDataset *AAIGDataset::CreateCopy(const char *pszFilename,
|
|
GDALDataset *poSrcDS, int /* bStrict */,
|
|
char **papszOptions,
|
|
GDALProgressFunc pfnProgress,
|
|
void *pProgressData)
|
|
{
|
|
const int nBands = poSrcDS->GetRasterCount();
|
|
const int nXSize = poSrcDS->GetRasterXSize();
|
|
const int nYSize = poSrcDS->GetRasterYSize();
|
|
|
|
// Some rudimentary checks.
|
|
if (nBands != 1)
|
|
{
|
|
ReportError(pszFilename, CE_Failure, CPLE_NotSupported,
|
|
"AAIG driver doesn't support %d bands. Must be 1 band.",
|
|
nBands);
|
|
|
|
return nullptr;
|
|
}
|
|
|
|
if (!pfnProgress(0.0, nullptr, pProgressData))
|
|
return nullptr;
|
|
|
|
// Create the dataset.
|
|
VSILFILE *fpImage = VSIFOpenL(pszFilename, "wt");
|
|
if (fpImage == nullptr)
|
|
{
|
|
ReportError(pszFilename, CE_Failure, CPLE_OpenFailed,
|
|
"Unable to create file.");
|
|
return nullptr;
|
|
}
|
|
|
|
// Write ASCII Grid file header.
|
|
double adfGeoTransform[6] = {};
|
|
char szHeader[2000] = {};
|
|
const char *pszForceCellsize =
|
|
CSLFetchNameValue(papszOptions, "FORCE_CELLSIZE");
|
|
|
|
poSrcDS->GetGeoTransform(adfGeoTransform);
|
|
|
|
const double dfYLLCorner =
|
|
adfGeoTransform[5] < 0
|
|
? adfGeoTransform[3] + nYSize * adfGeoTransform[5]
|
|
: adfGeoTransform[3];
|
|
if (std::abs(adfGeoTransform[1] + adfGeoTransform[5]) < 0.0000001 ||
|
|
std::abs(adfGeoTransform[1] - adfGeoTransform[5]) < 0.0000001 ||
|
|
(pszForceCellsize && CPLTestBool(pszForceCellsize)))
|
|
{
|
|
CPLsnprintf(szHeader, sizeof(szHeader),
|
|
"ncols %d\n"
|
|
"nrows %d\n"
|
|
"xllcorner %.12f\n"
|
|
"yllcorner %.12f\n"
|
|
"cellsize %.12f\n",
|
|
nXSize, nYSize, adfGeoTransform[0], dfYLLCorner,
|
|
adfGeoTransform[1]);
|
|
}
|
|
else
|
|
{
|
|
if (pszForceCellsize == nullptr)
|
|
ReportError(pszFilename, CE_Warning, CPLE_AppDefined,
|
|
"Producing a Golden Surfer style file with DX and DY "
|
|
"instead of CELLSIZE since the input pixels are "
|
|
"non-square. Use the FORCE_CELLSIZE=TRUE creation "
|
|
"option to force use of DX for even though this will "
|
|
"be distorted. Most ASCII Grid readers (ArcGIS "
|
|
"included) do not support the DX and DY parameters.");
|
|
CPLsnprintf(szHeader, sizeof(szHeader),
|
|
"ncols %d\n"
|
|
"nrows %d\n"
|
|
"xllcorner %.12f\n"
|
|
"yllcorner %.12f\n"
|
|
"dx %.12f\n"
|
|
"dy %.12f\n",
|
|
nXSize, nYSize, adfGeoTransform[0], dfYLLCorner,
|
|
adfGeoTransform[1], fabs(adfGeoTransform[5]));
|
|
}
|
|
|
|
// Builds the format string used for printing float values.
|
|
char szFormatFloat[32] = {'\0'};
|
|
strcpy(szFormatFloat, " %.20g");
|
|
const char *pszDecimalPrecision =
|
|
CSLFetchNameValue(papszOptions, "DECIMAL_PRECISION");
|
|
const char *pszSignificantDigits =
|
|
CSLFetchNameValue(papszOptions, "SIGNIFICANT_DIGITS");
|
|
bool bIgnoreSigDigits = false;
|
|
if (pszDecimalPrecision && pszSignificantDigits)
|
|
{
|
|
ReportError(pszFilename, CE_Warning, CPLE_AppDefined,
|
|
"Conflicting precision arguments, using DECIMAL_PRECISION");
|
|
bIgnoreSigDigits = true;
|
|
}
|
|
int nPrecision;
|
|
if (pszSignificantDigits && !bIgnoreSigDigits)
|
|
{
|
|
nPrecision = atoi(pszSignificantDigits);
|
|
if (nPrecision >= 0)
|
|
snprintf(szFormatFloat, sizeof(szFormatFloat), " %%.%dg",
|
|
nPrecision);
|
|
CPLDebug("AAIGrid", "Setting precision format: %s", szFormatFloat);
|
|
}
|
|
else if (pszDecimalPrecision)
|
|
{
|
|
nPrecision = atoi(pszDecimalPrecision);
|
|
if (nPrecision >= 0)
|
|
snprintf(szFormatFloat, sizeof(szFormatFloat), " %%.%df",
|
|
nPrecision);
|
|
CPLDebug("AAIGrid", "Setting precision format: %s", szFormatFloat);
|
|
}
|
|
|
|
// Handle nodata (optionally).
|
|
GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
|
|
const bool bReadAsInt = poBand->GetRasterDataType() == GDT_Byte ||
|
|
poBand->GetRasterDataType() == GDT_Int16 ||
|
|
poBand->GetRasterDataType() == GDT_UInt16 ||
|
|
poBand->GetRasterDataType() == GDT_Int32;
|
|
|
|
// Write `nodata' value to header if it is exists in source dataset
|
|
int bSuccess = FALSE;
|
|
const double dfNoData = poBand->GetNoDataValue(&bSuccess);
|
|
if (bSuccess)
|
|
{
|
|
snprintf(szHeader + strlen(szHeader),
|
|
sizeof(szHeader) - strlen(szHeader), "%s", "NODATA_value ");
|
|
if (bReadAsInt)
|
|
snprintf(szHeader + strlen(szHeader),
|
|
sizeof(szHeader) - strlen(szHeader), "%d",
|
|
static_cast<int>(dfNoData));
|
|
else
|
|
CPLsnprintf(szHeader + strlen(szHeader),
|
|
sizeof(szHeader) - strlen(szHeader), szFormatFloat,
|
|
dfNoData);
|
|
snprintf(szHeader + strlen(szHeader),
|
|
sizeof(szHeader) - strlen(szHeader), "%s", "\n");
|
|
}
|
|
|
|
if (VSIFWriteL(szHeader, strlen(szHeader), 1, fpImage) != 1)
|
|
{
|
|
CPL_IGNORE_RET_VAL(VSIFCloseL(fpImage));
|
|
return nullptr;
|
|
}
|
|
|
|
// Loop over image, copying image data.
|
|
|
|
// Write scanlines to output file
|
|
int *panScanline = bReadAsInt
|
|
? static_cast<int *>(CPLMalloc(
|
|
nXSize * GDALGetDataTypeSizeBytes(GDT_Int32)))
|
|
: nullptr;
|
|
|
|
double *padfScanline =
|
|
bReadAsInt ? nullptr
|
|
: static_cast<double *>(CPLMalloc(
|
|
nXSize * GDALGetDataTypeSizeBytes(GDT_Float64)));
|
|
|
|
CPLErr eErr = CE_None;
|
|
|
|
bool bHasOutputDecimalDot = false;
|
|
for (int iLine = 0; eErr == CE_None && iLine < nYSize; iLine++)
|
|
{
|
|
CPLString osBuf;
|
|
const int iSrcLine =
|
|
adfGeoTransform[5] < 0 ? iLine : nYSize - 1 - iLine;
|
|
eErr = poBand->RasterIO(GF_Read, 0, iSrcLine, nXSize, 1,
|
|
bReadAsInt ? static_cast<void *>(panScanline)
|
|
: static_cast<void *>(padfScanline),
|
|
nXSize, 1, bReadAsInt ? GDT_Int32 : GDT_Float64,
|
|
0, 0, nullptr);
|
|
|
|
if (bReadAsInt)
|
|
{
|
|
for (int iPixel = 0; iPixel < nXSize; iPixel++)
|
|
{
|
|
snprintf(szHeader, sizeof(szHeader), " %d",
|
|
panScanline[iPixel]);
|
|
osBuf += szHeader;
|
|
if ((iPixel & 1023) == 0 || iPixel == nXSize - 1)
|
|
{
|
|
if (VSIFWriteL(osBuf, static_cast<int>(osBuf.size()), 1,
|
|
fpImage) != 1)
|
|
{
|
|
eErr = CE_Failure;
|
|
ReportError(pszFilename, CE_Failure, CPLE_AppDefined,
|
|
"Write failed, disk full?");
|
|
break;
|
|
}
|
|
osBuf = "";
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
assert(padfScanline);
|
|
|
|
for (int iPixel = 0; iPixel < nXSize; iPixel++)
|
|
{
|
|
CPLsnprintf(szHeader, sizeof(szHeader), szFormatFloat,
|
|
padfScanline[iPixel]);
|
|
|
|
// Make sure that as least one value has a decimal point (#6060)
|
|
if (!bHasOutputDecimalDot)
|
|
{
|
|
if (strchr(szHeader, '.') || strchr(szHeader, 'e') ||
|
|
strchr(szHeader, 'E'))
|
|
{
|
|
bHasOutputDecimalDot = true;
|
|
}
|
|
else if (!CPLIsInf(padfScanline[iPixel]) &&
|
|
!CPLIsNan(padfScanline[iPixel]))
|
|
{
|
|
strcat(szHeader, ".0");
|
|
bHasOutputDecimalDot = true;
|
|
}
|
|
}
|
|
|
|
osBuf += szHeader;
|
|
if ((iPixel & 1023) == 0 || iPixel == nXSize - 1)
|
|
{
|
|
if (VSIFWriteL(osBuf, static_cast<int>(osBuf.size()), 1,
|
|
fpImage) != 1)
|
|
{
|
|
eErr = CE_Failure;
|
|
ReportError(pszFilename, CE_Failure, CPLE_AppDefined,
|
|
"Write failed, disk full?");
|
|
break;
|
|
}
|
|
osBuf = "";
|
|
}
|
|
}
|
|
}
|
|
if (VSIFWriteL("\n", 1, 1, fpImage) != 1)
|
|
eErr = CE_Failure;
|
|
|
|
if (eErr == CE_None &&
|
|
!pfnProgress((iLine + 1) / static_cast<double>(nYSize), nullptr,
|
|
pProgressData))
|
|
{
|
|
eErr = CE_Failure;
|
|
ReportError(pszFilename, CE_Failure, CPLE_UserInterrupt,
|
|
"User terminated CreateCopy()");
|
|
}
|
|
}
|
|
|
|
CPLFree(panScanline);
|
|
CPLFree(padfScanline);
|
|
if (VSIFCloseL(fpImage) != 0)
|
|
eErr = CE_Failure;
|
|
|
|
if (eErr != CE_None)
|
|
return nullptr;
|
|
|
|
// Try to write projection file.
|
|
const char *pszOriginalProjection = poSrcDS->GetProjectionRef();
|
|
if (!EQUAL(pszOriginalProjection, ""))
|
|
{
|
|
char *pszDirname = CPLStrdup(CPLGetPath(pszFilename));
|
|
char *pszBasename = CPLStrdup(CPLGetBasename(pszFilename));
|
|
char *pszPrjFilename =
|
|
CPLStrdup(CPLFormFilename(pszDirname, pszBasename, "prj"));
|
|
VSILFILE *fp = VSIFOpenL(pszPrjFilename, "wt");
|
|
if (fp != nullptr)
|
|
{
|
|
OGRSpatialReference oSRS;
|
|
oSRS.importFromWkt(pszOriginalProjection);
|
|
oSRS.morphToESRI();
|
|
char *pszESRIProjection = nullptr;
|
|
oSRS.exportToWkt(&pszESRIProjection);
|
|
CPL_IGNORE_RET_VAL(VSIFWriteL(pszESRIProjection, 1,
|
|
strlen(pszESRIProjection), fp));
|
|
|
|
CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
|
|
CPLFree(pszESRIProjection);
|
|
}
|
|
else
|
|
{
|
|
ReportError(pszFilename, CE_Failure, CPLE_FileIO,
|
|
"Unable to create file %s.", pszPrjFilename);
|
|
}
|
|
CPLFree(pszDirname);
|
|
CPLFree(pszBasename);
|
|
CPLFree(pszPrjFilename);
|
|
}
|
|
|
|
// Re-open dataset, and copy any auxiliary pam information.
|
|
|
|
// If writing to stdout, we can't reopen it, so return
|
|
// a fake dataset to make the caller happy.
|
|
CPLPushErrorHandler(CPLQuietErrorHandler);
|
|
GDALPamDataset *poDS =
|
|
reinterpret_cast<GDALPamDataset *>(GDALOpen(pszFilename, GA_ReadOnly));
|
|
CPLPopErrorHandler();
|
|
if (poDS)
|
|
{
|
|
poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
|
|
return poDS;
|
|
}
|
|
|
|
CPLErrorReset();
|
|
|
|
AAIGDataset *poAAIG_DS = new AAIGDataset();
|
|
poAAIG_DS->nRasterXSize = nXSize;
|
|
poAAIG_DS->nRasterYSize = nYSize;
|
|
poAAIG_DS->nBands = 1;
|
|
poAAIG_DS->SetBand(1, new AAIGRasterBand(poAAIG_DS, 1));
|
|
return poAAIG_DS;
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* OSR_GDS() */
|
|
/************************************************************************/
|
|
|
|
static CPLString OSR_GDS(char **papszNV, const char *pszField,
|
|
const char *pszDefaultValue)
|
|
|
|
{
|
|
if (papszNV == nullptr || papszNV[0] == nullptr)
|
|
return pszDefaultValue;
|
|
|
|
int iLine = 0; // Used after for.
|
|
for (; papszNV[iLine] != nullptr &&
|
|
!EQUALN(papszNV[iLine], pszField, strlen(pszField));
|
|
iLine++)
|
|
{
|
|
}
|
|
|
|
if (papszNV[iLine] == nullptr)
|
|
return pszDefaultValue;
|
|
else
|
|
{
|
|
char **papszTokens = CSLTokenizeString(papszNV[iLine]);
|
|
|
|
CPLString osResult;
|
|
if (CSLCount(papszTokens) > 1)
|
|
osResult = papszTokens[1];
|
|
else
|
|
osResult = pszDefaultValue;
|
|
|
|
CSLDestroy(papszTokens);
|
|
return osResult;
|
|
}
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GDALRegister_AAIGrid() */
|
|
/************************************************************************/
|
|
|
|
void GDALRegister_AAIGrid()
|
|
|
|
{
|
|
if (GDALGetDriverByName("AAIGrid") != nullptr)
|
|
return;
|
|
|
|
GDALDriver *poDriver = new GDALDriver();
|
|
|
|
poDriver->SetDescription("AAIGrid");
|
|
poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
|
|
poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Arc/Info ASCII Grid");
|
|
poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
|
|
"drivers/raster/aaigrid.html");
|
|
poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "asc");
|
|
poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
|
|
"Byte UInt16 Int16 Int32 Float32");
|
|
|
|
poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
|
|
poDriver->SetMetadataItem(
|
|
GDAL_DMD_CREATIONOPTIONLIST,
|
|
"<CreationOptionList>\n"
|
|
" <Option name='FORCE_CELLSIZE' type='boolean' description='Force "
|
|
"use of CELLSIZE, default is FALSE.'/>\n"
|
|
" <Option name='DECIMAL_PRECISION' type='int' description='Number of "
|
|
"decimal when writing floating-point numbers(%f).'/>\n"
|
|
" <Option name='SIGNIFICANT_DIGITS' type='int' description='Number "
|
|
"of significant digits when writing floating-point numbers(%g).'/>\n"
|
|
"</CreationOptionList>\n");
|
|
poDriver->SetMetadataItem(GDAL_DMD_OPENOPTIONLIST,
|
|
"<OpenOptionList>\n"
|
|
" <Option name='DATATYPE' type='string-select' "
|
|
"description='Data type to be used.'>\n"
|
|
" <Value>Int32</Value>\n"
|
|
" <Value>Float32</Value>\n"
|
|
" <Value>Float64</Value>\n"
|
|
" </Option>\n"
|
|
"</OpenOptionList>\n");
|
|
|
|
poDriver->pfnOpen = AAIGDataset::Open;
|
|
poDriver->pfnIdentify = AAIGDataset::Identify;
|
|
poDriver->pfnCreateCopy = AAIGDataset::CreateCopy;
|
|
|
|
GetGDALDriverManager()->RegisterDriver(poDriver);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GDALRegister_GRASSASCIIGrid() */
|
|
/************************************************************************/
|
|
|
|
void GDALRegister_GRASSASCIIGrid()
|
|
|
|
{
|
|
if (GDALGetDriverByName("GRASSASCIIGrid") != nullptr)
|
|
return;
|
|
|
|
GDALDriver *poDriver = new GDALDriver();
|
|
|
|
poDriver->SetDescription("GRASSASCIIGrid");
|
|
poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
|
|
poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "GRASS ASCII Grid");
|
|
poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
|
|
"drivers/raster/grassasciigrid.html");
|
|
|
|
poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
|
|
|
|
poDriver->pfnOpen = GRASSASCIIDataset::Open;
|
|
poDriver->pfnIdentify = GRASSASCIIDataset::Identify;
|
|
|
|
GetGDALDriverManager()->RegisterDriver(poDriver);
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* GDALRegister_ISG() */
|
|
/************************************************************************/
|
|
|
|
void GDALRegister_ISG()
|
|
|
|
{
|
|
if (GDALGetDriverByName("ISG") != nullptr)
|
|
return;
|
|
|
|
GDALDriver *poDriver = new GDALDriver();
|
|
|
|
poDriver->SetDescription("ISG");
|
|
poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
|
|
poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
|
|
"International Service for the Geoid");
|
|
poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/isg.html");
|
|
poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "isg");
|
|
|
|
poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
|
|
|
|
poDriver->pfnOpen = ISGDataset::Open;
|
|
poDriver->pfnIdentify = ISGDataset::Identify;
|
|
|
|
GetGDALDriverManager()->RegisterDriver(poDriver);
|
|
}
|