Coverage Report

Created: 2026-09-26 08:22

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/gdal/frmts/pcraster/pcrasterdataset.cpp
Line
Count
Source
1
/******************************************************************************
2
 *
3
 * Project:  PCRaster Integration
4
 * Purpose:  PCRaster CSF 2.0 raster file driver
5
 * Author:   Kor de Jong, Oliver Schmitz
6
 *
7
 ******************************************************************************
8
 * Copyright (c) PCRaster owners
9
 *
10
 * SPDX-License-Identifier: MIT
11
 ****************************************************************************/
12
13
#include "cpl_string.h"
14
#include "gdal_pam.h"
15
16
#include "pcrasterrasterband.h"
17
#include "pcrasterdataset.h"
18
#include "pcrasterutil.h"
19
#include "pcrasterdrivercore.h"
20
21
/*!
22
  \file
23
  This file contains the implementation of the PCRasterDataset class.
24
*/
25
26
//------------------------------------------------------------------------------
27
// DEFINITION OF STATIC PCRDATASET MEMBERS
28
//------------------------------------------------------------------------------
29
30
//! Tries to open the file described by \a info.
31
/*!
32
  \param     info Object with information about the dataset to open.
33
  \return    Pointer to newly allocated GDALDataset or 0.
34
35
  Returns a nullptr if the file could not be opened.
36
*/
37
GDALDataset *PCRasterDataset::open(GDALOpenInfo *info)
38
27
{
39
27
    PCRasterDataset *dataset = nullptr;
40
41
27
    if (PCRasterDriverIdentify(info))
42
27
    {
43
27
        MOPEN_PERM mode = info->eAccess == GA_Update ? M_READ_WRITE : M_READ;
44
45
27
        MAP *map = mapOpen(info->pszFilename, mode);
46
47
27
        if (map)
48
1
        {
49
1
            CPLErrorReset();
50
1
            dataset = new PCRasterDataset(map, info->eAccess);
51
1
            if (CPLGetLastErrorType() != CE_None)
52
0
            {
53
0
                delete dataset;
54
0
                return nullptr;
55
0
            }
56
1
        }
57
27
    }
58
59
    /* -------------------------------------------------------------------- */
60
    /*      Initialize any PAM information and overviews.                   */
61
    /* -------------------------------------------------------------------- */
62
27
    if (dataset)
63
1
    {
64
1
        dataset->SetDescription(info->pszFilename);
65
1
        dataset->TryLoadXML();
66
67
1
        dataset->oOvManager.Initialize(dataset, info->pszFilename);
68
1
    }
69
70
27
    return dataset;
71
27
}
72
73
//! Writes a raster to \a filename as a PCRaster raster file.
74
/*!
75
  \warning   The source raster must have only 1 band. Currently, the values in
76
             the source raster must be stored in one of the supported cell
77
             representations (CR_UINT1, CR_INT4, CR_REAL4, CR_REAL8).
78
79
  The meta data item PCRASTER_VALUESCALE will be checked to see what value
80
  scale to use. Otherwise a value scale is determined using
81
  GDALType2ValueScale(GDALDataType).
82
83
  This function always writes rasters using CR_UINT1, CR_INT4 or CR_REAL4
84
  cell representations.
85
*/
86
GDALDataset *PCRasterDataset::createCopy(char const *filename,
87
                                         GDALDataset *source,
88
                                         CPL_UNUSED int strict, CSLConstList,
89
                                         GDALProgressFunc progress,
90
                                         void *progressData)
91
0
{
92
    // Checks.
93
0
    const int nrBands = source->GetRasterCount();
94
0
    if (nrBands != 1)
95
0
    {
96
0
        CPLError(CE_Failure, CPLE_NotSupported,
97
0
                 "PCRaster driver: Too many bands ('%d'): must be 1 band",
98
0
                 nrBands);
99
0
        return nullptr;
100
0
    }
101
102
0
    GDALRasterBand *raster = source->GetRasterBand(1);
103
104
    // Create PCRaster raster. Determine properties of raster to create.
105
106
    // The in-file type of the cells.
107
0
    CSF_CR fileCellRepresentation =
108
0
        GDALType2CellRepresentation(raster->GetRasterDataType(), false);
109
110
0
    if (fileCellRepresentation == CR_UNDEFINED)
111
0
    {
112
0
        CPLError(
113
0
            CE_Failure, CPLE_NotSupported,
114
0
            "PCRaster driver: Cannot determine a valid cell representation");
115
0
        return nullptr;
116
0
    }
117
118
    // The value scale of the values.
119
0
    CSF_VS valueScale = VS_UNDEFINED;
120
0
    std::string osString;
121
0
    if (source->GetMetadataItem("PCRASTER_VALUESCALE"))
122
0
    {
123
0
        osString = source->GetMetadataItem("PCRASTER_VALUESCALE");
124
0
    }
125
126
0
    valueScale = !osString.empty()
127
0
                     ? string2ValueScale(osString)
128
0
                     : GDALType2ValueScale(raster->GetRasterDataType());
129
130
0
    if (valueScale == VS_UNDEFINED)
131
0
    {
132
0
        CPLError(CE_Failure, CPLE_NotSupported,
133
0
                 "PCRaster driver: Cannot determine a valid value scale");
134
0
        return nullptr;
135
0
    }
136
137
0
    CSF_PT const projection = PT_YDECT2B;
138
0
    const REAL8 angle = 0.0;
139
0
    REAL8 west = 0.0;
140
0
    REAL8 north = 0.0;
141
0
    REAL8 cellSize = 1.0;
142
143
0
    GDALGeoTransform gt;
144
0
    if (source->GetGeoTransform(gt) == CE_None)
145
0
    {
146
0
        if (gt.xrot == 0.0 && gt.yrot == 0.0)
147
0
        {
148
0
            west = static_cast<REAL8>(gt.xorig);
149
0
            north = static_cast<REAL8>(gt.yorig);
150
0
            cellSize = static_cast<REAL8>(gt.xscale);
151
0
        }
152
0
    }
153
154
    // The in-memory type of the cells.
155
0
    CSF_CR appCellRepresentation =
156
0
        GDALType2CellRepresentation(raster->GetRasterDataType(), true);
157
158
0
    if (appCellRepresentation == CR_UNDEFINED)
159
0
    {
160
0
        CPLError(
161
0
            CE_Failure, CPLE_NotSupported,
162
0
            "PCRaster driver: Cannot determine a valid cell representation");
163
0
        return nullptr;
164
0
    }
165
166
    // Check whether value scale fits the cell representation. Adjust when
167
    // needed.
168
0
    valueScale = fitValueScale(valueScale, appCellRepresentation);
169
170
    // Create a raster with the in file cell representation.
171
0
    const size_t nrRows = raster->GetYSize();
172
0
    const size_t nrCols = raster->GetXSize();
173
0
    MAP *map = Rcreate(filename, nrRows, nrCols, fileCellRepresentation,
174
0
                       valueScale, projection, west, north, angle, cellSize);
175
176
0
    if (!map)
177
0
    {
178
0
        CPLError(CE_Failure, CPLE_OpenFailed,
179
0
                 "PCRaster driver: Unable to create raster %s", filename);
180
0
        return nullptr;
181
0
    }
182
183
    // Try to convert in app cell representation to the cell representation
184
    // of the file.
185
0
    if (RuseAs(map, appCellRepresentation))
186
0
    {
187
0
        CPLError(CE_Failure, CPLE_NotSupported,
188
0
                 "PCRaster driver: Cannot convert cells: %s", MstrError());
189
0
        Mclose(map);
190
0
        return nullptr;
191
0
    }
192
193
0
    int hasMissingValue;
194
0
    double missingValue = raster->GetNoDataValue(&hasMissingValue);
195
196
    // This is needed to get my (KDJ) unit tests running.
197
    // I am still uncertain why this is needed. If the input raster has float32
198
    // values and the output int32, than the missing value in the dataset object
199
    // is not updated like the values are.
200
0
    if (missingValue == ::missingValue(CR_REAL4) &&
201
0
        fileCellRepresentation == CR_INT4)
202
0
    {
203
0
        missingValue = ::missingValue(fileCellRepresentation);
204
0
    }
205
206
    // TODO: Proper translation of TODO.
207
    // TODO: Support conversion to INT2 (?) INT4. ruseas.c see line 503.
208
    // Conversion r 159.
209
210
    // Create buffer for one row of values.
211
0
    void *buffer = Rmalloc(map, nrCols);
212
213
    // Copy values from source to target.
214
0
    CPLErr errorCode = CE_None;
215
0
    for (size_t row = 0; errorCode == CE_None && row < nrRows; ++row)
216
0
    {
217
218
        // Get row from source.
219
0
        if (raster->RasterIO(
220
0
                GF_Read, 0, static_cast<int>(row), static_cast<int>(nrCols), 1,
221
0
                buffer, static_cast<int>(nrCols), 1,
222
0
                raster->GetRasterDataType(), 0, 0, nullptr) != CE_None)
223
0
        {
224
0
            CPLError(CE_Failure, CPLE_FileIO,
225
0
                     "PCRaster driver: Error reading from source raster");
226
0
            errorCode = CE_Failure;
227
0
            break;
228
0
        }
229
230
        // Upon reading values are converted to the
231
        // right data type. This includes the missing value. If the source
232
        // value cannot be represented in the target data type it is set to a
233
        // missing value.
234
235
0
        if (hasMissingValue)
236
0
        {
237
0
            alterToStdMV(buffer, nrCols, appCellRepresentation, missingValue);
238
0
        }
239
240
0
        if (valueScale == VS_BOOLEAN)
241
0
        {
242
0
            castValuesToBooleanRange(buffer, nrCols, appCellRepresentation);
243
0
        }
244
245
        // Write row in target.
246
0
        RputRow(map, row, buffer);
247
248
0
        if (!progress((row + 1) / (static_cast<double>(nrRows)), nullptr,
249
0
                      progressData))
250
0
        {
251
0
            CPLError(CE_Failure, CPLE_UserInterrupt,
252
0
                     "PCRaster driver: User terminated CreateCopy()");
253
0
            errorCode = CE_Failure;
254
0
            break;
255
0
        }
256
0
    }
257
258
0
    Mclose(map);
259
0
    map = nullptr;
260
261
0
    free(buffer);
262
0
    buffer = nullptr;
263
264
0
    if (errorCode != CE_None)
265
0
        return nullptr;
266
267
    /* -------------------------------------------------------------------- */
268
    /*      Re-open dataset, and copy any auxiliary pam information.        */
269
    /* -------------------------------------------------------------------- */
270
0
    GDALPamDataset *poDS = cpl::down_cast<GDALPamDataset *>(
271
0
        GDALDataset::Open(filename, GDAL_OF_RASTER | GDAL_OF_UPDATE));
272
273
0
    if (poDS)
274
0
        poDS->CloneInfo(source, GCIF_PAM_DEFAULT);
275
276
0
    return poDS;
277
0
}
278
279
//------------------------------------------------------------------------------
280
// DEFINITION OF PCRDATASET MEMBERS
281
//------------------------------------------------------------------------------
282
283
//! Constructor.
284
/*!
285
  \param     mapIn PCRaster map handle. It is ours to close.
286
*/
287
PCRasterDataset::PCRasterDataset(MAP *mapIn, GDALAccess eAccessIn)
288
1
    : GDALPamDataset(), d_map(mapIn), d_west(0.0), d_north(0.0),
289
1
      d_cellSize(0.0), d_cellRepresentation(CR_UNDEFINED),
290
1
      d_valueScale(VS_UNDEFINED), d_defaultNoDataValue(0.0),
291
1
      d_location_changed(false)
292
1
{
293
    // Read header info.
294
1
    eAccess = eAccessIn;
295
1
    nRasterXSize = static_cast<int>(RgetNrCols(d_map));
296
1
    nRasterYSize = static_cast<int>(RgetNrRows(d_map));
297
1
    if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
298
0
    {
299
0
        return;
300
0
    }
301
1
    d_west = static_cast<double>(RgetXUL(d_map));
302
1
    d_north = static_cast<double>(RgetYUL(d_map));
303
1
    d_cellSize = static_cast<double>(RgetCellSize(d_map));
304
1
    d_cellRepresentation = RgetUseCellRepr(d_map);
305
1
    if (d_cellRepresentation == CR_UNDEFINED)
306
0
    {
307
0
        CPLError(CE_Failure, CPLE_AssertionFailed,
308
0
                 "d_cellRepresentation != CR_UNDEFINED");
309
0
    }
310
1
    d_valueScale = RgetValueScale(d_map);
311
1
    if (d_valueScale == VS_UNDEFINED)
312
0
    {
313
0
        CPLError(CE_Failure, CPLE_AssertionFailed,
314
0
                 "d_valueScale != VS_UNDEFINED");
315
0
    }
316
1
    d_defaultNoDataValue = ::missingValue(d_cellRepresentation);
317
318
    // Create band information objects.
319
1
    nBands = 1;
320
1
    SetBand(1, new PCRasterRasterBand(this));
321
322
1
    SetMetadataItem("PCRASTER_VALUESCALE",
323
1
                    valueScale2String(d_valueScale).c_str());
324
1
}
325
326
//! Destructor.
327
/*!
328
  \warning   The map given in the constructor is closed.
329
*/
330
PCRasterDataset::~PCRasterDataset()
331
1
{
332
1
    FlushCache(true);
333
1
    Mclose(d_map);
334
1
}
335
336
//! Sets projections info.
337
/*!
338
  \param     gt Array to fill.
339
340
  CSF 2.0 supports the notion of y coordinates which increase from north to
341
  south. Support for this has been dropped and applications reading PCRaster
342
  rasters will treat or already treat y coordinates as increasing from south
343
  to north only.
344
*/
345
CPLErr PCRasterDataset::GetGeoTransform(GDALGeoTransform &gt) const
346
1
{
347
    // x = west + nrCols * cellsize
348
1
    gt.xorig = d_west;
349
1
    gt.xscale = d_cellSize;
350
1
    gt.xrot = 0.0;
351
352
    // y = north + nrRows * -cellsize
353
1
    gt.yorig = d_north;
354
1
    gt.yrot = 0.0;
355
1
    gt.yscale = -1.0 * d_cellSize;
356
357
1
    return CE_None;
358
1
}
359
360
//! Returns the map handle.
361
/*!
362
  \return    Map handle.
363
*/
364
MAP *PCRasterDataset::map() const
365
512
{
366
512
    return d_map;
367
512
}
368
369
//! Returns the in-app cell representation.
370
/*!
371
  \return    cell representation
372
  \warning   This might not be the same representation as use to store the
373
  values in the file. \sa        valueScale()
374
*/
375
CSF_CR PCRasterDataset::cellRepresentation() const
376
513
{
377
513
    return d_cellRepresentation;
378
513
}
379
380
//! Returns the value scale of the data.
381
/*!
382
  \return    Value scale
383
  \sa        cellRepresentation()
384
*/
385
CSF_VS PCRasterDataset::valueScale() const
386
0
{
387
0
    return d_valueScale;
388
0
}
389
390
//! Returns the value of the missing value.
391
/*!
392
  \return    Missing value
393
*/
394
double PCRasterDataset::defaultNoDataValue() const
395
516
{
396
516
    return d_defaultNoDataValue;
397
516
}
398
399
GDALDataset *PCRasterDataset::create(const char *filename, int nr_cols,
400
                                     int nr_rows, int nrBands,
401
                                     GDALDataType gdalType,
402
                                     CSLConstList papszParamList)
403
0
{
404
    // Checks
405
0
    if (nrBands != 1)
406
0
    {
407
0
        CPLError(CE_Failure, CPLE_NotSupported,
408
0
                 "PCRaster driver : "
409
0
                 "attempt to create dataset with too many bands (%d); "
410
0
                 "must be 1 band.\n",
411
0
                 nrBands);
412
0
        return nullptr;
413
0
    }
414
415
0
    const int row_col_max = INT4_MAX - 1;
416
0
    if (nr_cols > row_col_max)
417
0
    {
418
0
        CPLError(CE_Failure, CPLE_NotSupported,
419
0
                 "PCRaster driver : "
420
0
                 "attempt to create dataset with too many columns (%d); "
421
0
                 "must be smaller than %d.",
422
0
                 nr_cols, row_col_max);
423
0
        return nullptr;
424
0
    }
425
426
0
    if (nr_rows > row_col_max)
427
0
    {
428
0
        CPLError(CE_Failure, CPLE_NotSupported,
429
0
                 "PCRaster driver : "
430
0
                 "attempt to create dataset with too many rows (%d); "
431
0
                 "must be smaller than %d.",
432
0
                 nr_rows, row_col_max);
433
0
        return nullptr;
434
0
    }
435
436
0
    if (gdalType != GDT_UInt8 && gdalType != GDT_Int32 &&
437
0
        gdalType != GDT_Float32)
438
0
    {
439
0
        CPLError(CE_Failure, CPLE_AppDefined,
440
0
                 "PCRaster driver: "
441
0
                 "attempt to create dataset with an illegal data type (%s); "
442
0
                 "use either Byte, Int32 or Float32.",
443
0
                 GDALGetDataTypeName(gdalType));
444
0
        return nullptr;
445
0
    }
446
447
    // value scale must be specified by the user,
448
    // determines cell representation
449
0
    const char *valueScale =
450
0
        CSLFetchNameValue(papszParamList, "PCRASTER_VALUESCALE");
451
452
0
    if (valueScale == nullptr)
453
0
    {
454
0
        CPLError(CE_Failure, CPLE_AppDefined,
455
0
                 "PCRaster driver: value scale can not be determined; "
456
0
                 "specify PCRASTER_VALUESCALE.");
457
0
        return nullptr;
458
0
    }
459
460
0
    CSF_VS csf_value_scale = string2ValueScale(valueScale);
461
462
0
    if (csf_value_scale == VS_UNDEFINED)
463
0
    {
464
0
        CPLError(CE_Failure, CPLE_AppDefined,
465
0
                 "PCRaster driver: value scale can not be determined (%s); "
466
0
                 "use either VS_BOOLEAN, VS_NOMINAL, VS_ORDINAL, VS_SCALAR, "
467
0
                 "VS_DIRECTION, VS_LDD",
468
0
                 valueScale);
469
0
        return nullptr;
470
0
    }
471
472
0
    CSF_CR csf_cell_representation =
473
0
        GDALType2CellRepresentation(gdalType, false);
474
475
    // default values
476
0
    REAL8 west = 0.0;
477
0
    REAL8 north = 0.0;
478
0
    REAL8 length = 1.0;
479
0
    REAL8 angle = 0.0;
480
0
    CSF_PT projection = PT_YDECT2B;
481
482
    // Create a new raster
483
0
    MAP *map = Rcreate(filename, nr_rows, nr_cols, csf_cell_representation,
484
0
                       csf_value_scale, projection, west, north, angle, length);
485
486
0
    if (!map)
487
0
    {
488
0
        CPLError(CE_Failure, CPLE_OpenFailed,
489
0
                 "PCRaster driver: Unable to create raster %s", filename);
490
0
        return nullptr;
491
0
    }
492
493
0
    Mclose(map);
494
0
    map = nullptr;
495
496
0
    return GDALDataset::Open(filename, GDAL_OF_RASTER | GDAL_OF_UPDATE);
497
0
}
498
499
CPLErr PCRasterDataset::SetGeoTransform(const GDALGeoTransform &gt)
500
0
{
501
0
    if ((gt.xrot != 0.0) || (gt.yrot != 0.0))
502
0
    {
503
0
        CPLError(CE_Failure, CPLE_NotSupported,
504
0
                 "PCRaster driver: "
505
0
                 "rotated geogtations are not supported.");
506
0
        return CE_Failure;
507
0
    }
508
509
0
    if (gt.xscale != gt.yscale * -1.0)
510
0
    {
511
0
        CPLError(CE_Failure, CPLE_NotSupported,
512
0
                 "PCRaster driver: "
513
0
                 "only the same width and height for cells is supported.");
514
0
        return CE_Failure;
515
0
    }
516
517
0
    d_west = gt.xorig;
518
0
    d_north = gt.yorig;
519
0
    d_cellSize = gt.xscale;
520
0
    d_location_changed = true;
521
522
0
    return CE_None;
523
0
}
524
525
bool PCRasterDataset::location_changed() const
526
0
{
527
0
    return d_location_changed;
528
0
}