Coverage Report

Created: 2025-06-09 07:07

/src/gdal/frmts/jdem/jdemdataset.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  JDEM Reader
4
 * Purpose:  All code for Japanese DEM Reader
5
 * Author:   Frank Warmerdam, warmerdam@pobox.com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
9
 * Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 ****************************************************************************/
13
14
#include "cpl_port.h"
15
#include "gdal_frmts.h"
16
#include "gdal_pam.h"
17
18
#include <algorithm>
19
20
constexpr int HEADER_SIZE = 1011;
21
22
/************************************************************************/
23
/*                            JDEMGetField()                            */
24
/************************************************************************/
25
26
static int JDEMGetField(const char *pszField, int nWidth)
27
28
0
{
29
0
    char szWork[32] = {};
30
0
    CPLAssert(nWidth < static_cast<int>(sizeof(szWork)));
31
32
0
    strncpy(szWork, pszField, nWidth);
33
0
    szWork[nWidth] = '\0';
34
35
0
    return atoi(szWork);
36
0
}
37
38
/************************************************************************/
39
/*                            JDEMGetAngle()                            */
40
/************************************************************************/
41
42
static double JDEMGetAngle(const char *pszField)
43
44
0
{
45
0
    const int nAngle = JDEMGetField(pszField, 7);
46
47
    // Note, this isn't very general purpose, but it would appear
48
    // from the field widths that angles are never negative.  Nice
49
    // to be a country in the "first quadrant".
50
51
0
    const int nDegree = nAngle / 10000;
52
0
    const int nMin = (nAngle / 100) % 100;
53
0
    const int nSec = nAngle % 100;
54
55
0
    return nDegree + nMin / 60.0 + nSec / 3600.0;
56
0
}
57
58
/************************************************************************/
59
/* ==================================================================== */
60
/*                              JDEMDataset                             */
61
/* ==================================================================== */
62
/************************************************************************/
63
64
class JDEMRasterBand;
65
66
class JDEMDataset final : public GDALPamDataset
67
{
68
    friend class JDEMRasterBand;
69
70
    VSILFILE *m_fp = nullptr;
71
    GByte m_abyHeader[HEADER_SIZE];
72
    OGRSpatialReference m_oSRS{};
73
74
  public:
75
    JDEMDataset();
76
    ~JDEMDataset();
77
78
    static GDALDataset *Open(GDALOpenInfo *);
79
    static int Identify(GDALOpenInfo *);
80
81
    CPLErr GetGeoTransform(double *padfTransform) override;
82
    const OGRSpatialReference *GetSpatialRef() const override;
83
};
84
85
/************************************************************************/
86
/* ==================================================================== */
87
/*                            JDEMRasterBand                             */
88
/* ==================================================================== */
89
/************************************************************************/
90
91
class JDEMRasterBand final : public GDALPamRasterBand
92
{
93
    friend class JDEMDataset;
94
95
    int m_nRecordSize = 0;
96
    char *m_pszRecord = nullptr;
97
    bool m_bBufferAllocFailed = false;
98
99
  public:
100
    JDEMRasterBand(JDEMDataset *, int);
101
    ~JDEMRasterBand();
102
103
    virtual CPLErr IReadBlock(int, int, void *) override;
104
};
105
106
/************************************************************************/
107
/*                           JDEMRasterBand()                            */
108
/************************************************************************/
109
110
JDEMRasterBand::JDEMRasterBand(JDEMDataset *poDSIn, int nBandIn)
111
    :  // Cannot overflow as nBlockXSize <= 999.
112
0
      m_nRecordSize(poDSIn->GetRasterXSize() * 5 + 9 + 2)
113
0
{
114
0
    poDS = poDSIn;
115
0
    nBand = nBandIn;
116
117
0
    eDataType = GDT_Float32;
118
119
0
    nBlockXSize = poDS->GetRasterXSize();
120
0
    nBlockYSize = 1;
121
0
}
122
123
/************************************************************************/
124
/*                          ~JDEMRasterBand()                            */
125
/************************************************************************/
126
127
JDEMRasterBand::~JDEMRasterBand()
128
0
{
129
0
    VSIFree(m_pszRecord);
130
0
}
131
132
/************************************************************************/
133
/*                             IReadBlock()                             */
134
/************************************************************************/
135
136
CPLErr JDEMRasterBand::IReadBlock(int /* nBlockXOff */, int nBlockYOff,
137
                                  void *pImage)
138
139
0
{
140
0
    JDEMDataset *poGDS = cpl::down_cast<JDEMDataset *>(poDS);
141
142
0
    if (m_pszRecord == nullptr)
143
0
    {
144
0
        if (m_bBufferAllocFailed)
145
0
            return CE_Failure;
146
147
0
        m_pszRecord = static_cast<char *>(VSI_MALLOC_VERBOSE(m_nRecordSize));
148
0
        if (m_pszRecord == nullptr)
149
0
        {
150
0
            m_bBufferAllocFailed = true;
151
0
            return CE_Failure;
152
0
        }
153
0
    }
154
155
0
    CPL_IGNORE_RET_VAL(
156
0
        VSIFSeekL(poGDS->m_fp, 1011 + m_nRecordSize * nBlockYOff, SEEK_SET));
157
158
0
    if (VSIFReadL(m_pszRecord, m_nRecordSize, 1, poGDS->m_fp) != 1)
159
0
    {
160
0
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot read scanline %d",
161
0
                 nBlockYOff);
162
0
        return CE_Failure;
163
0
    }
164
165
0
    if (!EQUALN(reinterpret_cast<char *>(poGDS->m_abyHeader), m_pszRecord, 6))
166
0
    {
167
0
        CPLError(CE_Failure, CPLE_AppDefined,
168
0
                 "JDEM Scanline corrupt.  Perhaps file was not transferred "
169
0
                 "in binary mode?");
170
0
        return CE_Failure;
171
0
    }
172
173
0
    if (JDEMGetField(m_pszRecord + 6, 3) != nBlockYOff + 1)
174
0
    {
175
0
        CPLError(CE_Failure, CPLE_AppDefined,
176
0
                 "JDEM scanline out of order, JDEM driver does not "
177
0
                 "currently support partial datasets.");
178
0
        return CE_Failure;
179
0
    }
180
181
0
    for (int i = 0; i < nBlockXSize; i++)
182
0
        static_cast<float *>(pImage)[i] =
183
0
            JDEMGetField(m_pszRecord + 9 + 5 * i, 5) * 0.1f;
184
185
0
    return CE_None;
186
0
}
187
188
/************************************************************************/
189
/* ==================================================================== */
190
/*                              JDEMDataset                             */
191
/* ==================================================================== */
192
/************************************************************************/
193
194
/************************************************************************/
195
/*                            JDEMDataset()                             */
196
/************************************************************************/
197
198
JDEMDataset::JDEMDataset()
199
0
{
200
0
    std::fill_n(m_abyHeader, CPL_ARRAYSIZE(m_abyHeader), static_cast<GByte>(0));
201
0
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
202
0
    m_oSRS.importFromEPSG(4301);  // Tokyo geographic CRS
203
0
}
204
205
/************************************************************************/
206
/*                           ~JDEMDataset()                             */
207
/************************************************************************/
208
209
JDEMDataset::~JDEMDataset()
210
211
0
{
212
0
    FlushCache(true);
213
0
    if (m_fp != nullptr)
214
0
        CPL_IGNORE_RET_VAL(VSIFCloseL(m_fp));
215
0
}
216
217
/************************************************************************/
218
/*                          GetGeoTransform()                           */
219
/************************************************************************/
220
221
CPLErr JDEMDataset::GetGeoTransform(double *padfTransform)
222
223
0
{
224
0
    const char *psHeader = reinterpret_cast<const char *>(m_abyHeader);
225
226
0
    const double dfLLLat = JDEMGetAngle(psHeader + 29);
227
0
    const double dfLLLong = JDEMGetAngle(psHeader + 36);
228
0
    const double dfURLat = JDEMGetAngle(psHeader + 43);
229
0
    const double dfURLong = JDEMGetAngle(psHeader + 50);
230
231
0
    padfTransform[0] = dfLLLong;
232
0
    padfTransform[3] = dfURLat;
233
0
    padfTransform[1] = (dfURLong - dfLLLong) / GetRasterXSize();
234
0
    padfTransform[2] = 0.0;
235
236
0
    padfTransform[4] = 0.0;
237
0
    padfTransform[5] = -1 * (dfURLat - dfLLLat) / GetRasterYSize();
238
239
0
    return CE_None;
240
0
}
241
242
/************************************************************************/
243
/*                          GetSpatialRef()                             */
244
/************************************************************************/
245
246
const OGRSpatialReference *JDEMDataset::GetSpatialRef() const
247
248
0
{
249
0
    return &m_oSRS;
250
0
}
251
252
/************************************************************************/
253
/*                              Identify()                              */
254
/************************************************************************/
255
256
int JDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
257
258
16.3k
{
259
16.3k
    if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
260
3.95k
        return FALSE;
261
262
    // Confirm that the header has what appears to be dates in the
263
    // expected locations.
264
    // Check if century values seem reasonable.
265
12.3k
    const char *psHeader = reinterpret_cast<char *>(poOpenInfo->pabyHeader);
266
12.3k
    if ((!STARTS_WITH_CI(psHeader + 11, "19") &&
267
12.3k
         !STARTS_WITH_CI(psHeader + 11, "20")) ||
268
12.3k
        (!STARTS_WITH_CI(psHeader + 15, "19") &&
269
28
         !STARTS_WITH_CI(psHeader + 15, "20")) ||
270
12.3k
        (!STARTS_WITH_CI(psHeader + 19, "19") &&
271
0
         !STARTS_WITH_CI(psHeader + 19, "20")))
272
12.3k
    {
273
12.3k
        return FALSE;
274
12.3k
    }
275
276
    // Check the extent too. In particular, that we are in the first quadrant,
277
    // as this is only for Japan.
278
0
    const double dfLLLat = JDEMGetAngle(psHeader + 29);
279
0
    const double dfLLLong = JDEMGetAngle(psHeader + 36);
280
0
    const double dfURLat = JDEMGetAngle(psHeader + 43);
281
0
    const double dfURLong = JDEMGetAngle(psHeader + 50);
282
0
    if (dfLLLat > 90 || dfLLLat < 0 || dfLLLong > 180 || dfLLLong < 0 ||
283
0
        dfURLat > 90 || dfURLat < 0 || dfURLong > 180 || dfURLong < 0 ||
284
0
        dfLLLat > dfURLat || dfLLLong > dfURLong)
285
0
    {
286
0
        return FALSE;
287
0
    }
288
289
0
    return TRUE;
290
0
}
291
292
/************************************************************************/
293
/*                                Open()                                */
294
/************************************************************************/
295
296
GDALDataset *JDEMDataset::Open(GDALOpenInfo *poOpenInfo)
297
298
0
{
299
    // Confirm that the header is compatible with a JDEM dataset.
300
0
    if (!Identify(poOpenInfo))
301
0
        return nullptr;
302
303
    // Confirm the requested access is supported.
304
0
    if (poOpenInfo->eAccess == GA_Update)
305
0
    {
306
0
        ReportUpdateNotSupportedByDriver("JDEM");
307
0
        return nullptr;
308
0
    }
309
310
    // Check that the file pointer from GDALOpenInfo* is available.
311
0
    if (poOpenInfo->fpL == nullptr)
312
0
    {
313
0
        return nullptr;
314
0
    }
315
316
    // Create a corresponding GDALDataset.
317
0
    auto poDS = std::make_unique<JDEMDataset>();
318
319
    // Borrow the file pointer from GDALOpenInfo*.
320
0
    std::swap(poDS->m_fp, poOpenInfo->fpL);
321
322
    // Store the header (we have already checked it is at least HEADER_SIZE
323
    // byte large).
324
0
    memcpy(poDS->m_abyHeader, poOpenInfo->pabyHeader, HEADER_SIZE);
325
326
0
    const char *psHeader = reinterpret_cast<const char *>(poDS->m_abyHeader);
327
0
    poDS->nRasterXSize = JDEMGetField(psHeader + 23, 3);
328
0
    poDS->nRasterYSize = JDEMGetField(psHeader + 26, 3);
329
0
    if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
330
0
    {
331
0
        return nullptr;
332
0
    }
333
334
    // Create band information objects.
335
0
    poDS->SetBand(1, new JDEMRasterBand(poDS.get(), 1));
336
337
    // Initialize any PAM information.
338
0
    poDS->SetDescription(poOpenInfo->pszFilename);
339
0
    poDS->TryLoadXML();
340
341
    // Check for overviews.
342
0
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
343
344
0
    return poDS.release();
345
0
}
346
347
/************************************************************************/
348
/*                          GDALRegister_JDEM()                         */
349
/************************************************************************/
350
351
void GDALRegister_JDEM()
352
353
2
{
354
2
    if (GDALGetDriverByName("JDEM") != nullptr)
355
0
        return;
356
357
2
    GDALDriver *poDriver = new GDALDriver();
358
359
2
    poDriver->SetDescription("JDEM");
360
2
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
361
2
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Japanese DEM (.mem)");
362
2
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/jdem.html");
363
2
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "mem");
364
2
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
365
366
2
    poDriver->pfnOpen = JDEMDataset::Open;
367
2
    poDriver->pfnIdentify = JDEMDataset::Identify;
368
369
2
    GetGDALDriverManager()->RegisterDriver(poDriver);
370
2
}