Line data Source code
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 13 : GDALDataset *PCRasterDataset::open(GDALOpenInfo *info)
38 : {
39 13 : PCRasterDataset *dataset = nullptr;
40 :
41 13 : if (PCRasterDriverIdentify(info))
42 : {
43 13 : MOPEN_PERM mode = info->eAccess == GA_Update ? M_READ_WRITE : M_READ;
44 :
45 13 : MAP *map = mapOpen(info->pszFilename, mode);
46 :
47 13 : if (map)
48 : {
49 13 : CPLErrorReset();
50 13 : dataset = new PCRasterDataset(map, info->eAccess);
51 13 : if (CPLGetLastErrorType() != CE_None)
52 : {
53 0 : delete dataset;
54 0 : return nullptr;
55 : }
56 : }
57 : }
58 :
59 : /* -------------------------------------------------------------------- */
60 : /* Initialize any PAM information and overviews. */
61 : /* -------------------------------------------------------------------- */
62 13 : if (dataset)
63 : {
64 13 : dataset->SetDescription(info->pszFilename);
65 13 : dataset->TryLoadXML();
66 :
67 13 : dataset->oOvManager.Initialize(dataset, info->pszFilename);
68 : }
69 :
70 13 : return dataset;
71 : }
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 20 : GDALDataset *PCRasterDataset::createCopy(char const *filename,
87 : GDALDataset *source,
88 : CPL_UNUSED int strict, CSLConstList,
89 : GDALProgressFunc progress,
90 : void *progressData)
91 : {
92 : // Checks.
93 20 : const int nrBands = source->GetRasterCount();
94 20 : if (nrBands != 1)
95 : {
96 5 : CPLError(CE_Failure, CPLE_NotSupported,
97 : "PCRaster driver: Too many bands ('%d'): must be 1 band",
98 : nrBands);
99 5 : return nullptr;
100 : }
101 :
102 15 : 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 : CSF_CR fileCellRepresentation =
108 15 : GDALType2CellRepresentation(raster->GetRasterDataType(), false);
109 :
110 15 : if (fileCellRepresentation == CR_UNDEFINED)
111 : {
112 4 : CPLError(
113 : CE_Failure, CPLE_NotSupported,
114 : "PCRaster driver: Cannot determine a valid cell representation");
115 4 : return nullptr;
116 : }
117 :
118 : // The value scale of the values.
119 11 : CSF_VS valueScale = VS_UNDEFINED;
120 22 : std::string osString;
121 11 : if (source->GetMetadataItem("PCRASTER_VALUESCALE"))
122 : {
123 1 : osString = source->GetMetadataItem("PCRASTER_VALUESCALE");
124 : }
125 :
126 11 : valueScale = !osString.empty()
127 11 : ? string2ValueScale(osString)
128 10 : : GDALType2ValueScale(raster->GetRasterDataType());
129 :
130 11 : if (valueScale == VS_UNDEFINED)
131 : {
132 0 : CPLError(CE_Failure, CPLE_NotSupported,
133 : "PCRaster driver: Cannot determine a valid value scale");
134 0 : return nullptr;
135 : }
136 :
137 11 : CSF_PT const projection = PT_YDECT2B;
138 11 : const REAL8 angle = 0.0;
139 11 : REAL8 west = 0.0;
140 11 : REAL8 north = 0.0;
141 11 : REAL8 cellSize = 1.0;
142 :
143 11 : GDALGeoTransform gt;
144 11 : if (source->GetGeoTransform(gt) == CE_None)
145 : {
146 11 : if (gt[2] == 0.0 && gt[4] == 0.0)
147 : {
148 11 : west = static_cast<REAL8>(gt[0]);
149 11 : north = static_cast<REAL8>(gt[3]);
150 11 : cellSize = static_cast<REAL8>(gt[1]);
151 : }
152 : }
153 :
154 : // The in-memory type of the cells.
155 : CSF_CR appCellRepresentation =
156 11 : GDALType2CellRepresentation(raster->GetRasterDataType(), true);
157 :
158 11 : if (appCellRepresentation == CR_UNDEFINED)
159 : {
160 0 : CPLError(
161 : CE_Failure, CPLE_NotSupported,
162 : "PCRaster driver: Cannot determine a valid cell representation");
163 0 : return nullptr;
164 : }
165 :
166 : // Check whether value scale fits the cell representation. Adjust when
167 : // needed.
168 11 : valueScale = fitValueScale(valueScale, appCellRepresentation);
169 :
170 : // Create a raster with the in file cell representation.
171 11 : const size_t nrRows = raster->GetYSize();
172 11 : const size_t nrCols = raster->GetXSize();
173 11 : MAP *map = Rcreate(filename, nrRows, nrCols, fileCellRepresentation,
174 : valueScale, projection, west, north, angle, cellSize);
175 :
176 11 : if (!map)
177 : {
178 3 : CPLError(CE_Failure, CPLE_OpenFailed,
179 : "PCRaster driver: Unable to create raster %s", filename);
180 3 : return nullptr;
181 : }
182 :
183 : // Try to convert in app cell representation to the cell representation
184 : // of the file.
185 8 : if (RuseAs(map, appCellRepresentation))
186 : {
187 3 : CPLError(CE_Failure, CPLE_NotSupported,
188 : "PCRaster driver: Cannot convert cells: %s", MstrError());
189 3 : Mclose(map);
190 3 : return nullptr;
191 : }
192 :
193 : int hasMissingValue;
194 5 : 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 5 : if (missingValue == ::missingValue(CR_REAL4) &&
201 : fileCellRepresentation == CR_INT4)
202 : {
203 0 : missingValue = ::missingValue(fileCellRepresentation);
204 : }
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 5 : void *buffer = Rmalloc(map, nrCols);
212 :
213 : // Copy values from source to target.
214 5 : CPLErr errorCode = CE_None;
215 145 : for (size_t row = 0; errorCode == CE_None && row < nrRows; ++row)
216 : {
217 :
218 : // Get row from source.
219 140 : if (raster->RasterIO(
220 : GF_Read, 0, static_cast<int>(row), static_cast<int>(nrCols), 1,
221 : buffer, static_cast<int>(nrCols), 1,
222 140 : raster->GetRasterDataType(), 0, 0, nullptr) != CE_None)
223 : {
224 0 : CPLError(CE_Failure, CPLE_FileIO,
225 : "PCRaster driver: Error reading from source raster");
226 0 : errorCode = CE_Failure;
227 0 : break;
228 : }
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 140 : if (hasMissingValue)
236 : {
237 100 : alterToStdMV(buffer, nrCols, appCellRepresentation, missingValue);
238 : }
239 :
240 140 : if (valueScale == VS_BOOLEAN)
241 : {
242 10 : castValuesToBooleanRange(buffer, nrCols, appCellRepresentation);
243 : }
244 :
245 : // Write row in target.
246 140 : RputRow(map, row, buffer);
247 :
248 140 : if (!progress((row + 1) / (static_cast<double>(nrRows)), nullptr,
249 : progressData))
250 : {
251 0 : CPLError(CE_Failure, CPLE_UserInterrupt,
252 : "PCRaster driver: User terminated CreateCopy()");
253 0 : errorCode = CE_Failure;
254 0 : break;
255 : }
256 : }
257 :
258 5 : Mclose(map);
259 5 : map = nullptr;
260 :
261 5 : free(buffer);
262 5 : buffer = nullptr;
263 :
264 5 : if (errorCode != CE_None)
265 0 : return nullptr;
266 :
267 : /* -------------------------------------------------------------------- */
268 : /* Re-open dataset, and copy any auxiliary pam information. */
269 : /* -------------------------------------------------------------------- */
270 5 : GDALPamDataset *poDS = cpl::down_cast<GDALPamDataset *>(
271 : GDALDataset::Open(filename, GDAL_OF_RASTER | GDAL_OF_UPDATE));
272 :
273 5 : if (poDS)
274 5 : poDS->CloneInfo(source, GCIF_PAM_DEFAULT);
275 :
276 5 : return poDS;
277 : }
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 13 : PCRasterDataset::PCRasterDataset(MAP *mapIn, GDALAccess eAccessIn)
288 : : GDALPamDataset(), d_map(mapIn), d_west(0.0), d_north(0.0),
289 : d_cellSize(0.0), d_cellRepresentation(CR_UNDEFINED),
290 : d_valueScale(VS_UNDEFINED), d_defaultNoDataValue(0.0),
291 13 : d_location_changed(false)
292 : {
293 : // Read header info.
294 13 : eAccess = eAccessIn;
295 13 : nRasterXSize = static_cast<int>(RgetNrCols(d_map));
296 13 : nRasterYSize = static_cast<int>(RgetNrRows(d_map));
297 13 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
298 : {
299 0 : return;
300 : }
301 13 : d_west = static_cast<double>(RgetXUL(d_map));
302 13 : d_north = static_cast<double>(RgetYUL(d_map));
303 13 : d_cellSize = static_cast<double>(RgetCellSize(d_map));
304 13 : d_cellRepresentation = RgetUseCellRepr(d_map);
305 13 : if (d_cellRepresentation == CR_UNDEFINED)
306 : {
307 0 : CPLError(CE_Failure, CPLE_AssertionFailed,
308 : "d_cellRepresentation != CR_UNDEFINED");
309 : }
310 13 : d_valueScale = RgetValueScale(d_map);
311 13 : if (d_valueScale == VS_UNDEFINED)
312 : {
313 0 : CPLError(CE_Failure, CPLE_AssertionFailed,
314 : "d_valueScale != VS_UNDEFINED");
315 : }
316 13 : d_defaultNoDataValue = ::missingValue(d_cellRepresentation);
317 :
318 : // Create band information objects.
319 13 : nBands = 1;
320 13 : SetBand(1, new PCRasterRasterBand(this));
321 :
322 13 : SetMetadataItem("PCRASTER_VALUESCALE",
323 26 : valueScale2String(d_valueScale).c_str());
324 : }
325 :
326 : //! Destructor.
327 : /*!
328 : \warning The map given in the constructor is closed.
329 : */
330 26 : PCRasterDataset::~PCRasterDataset()
331 : {
332 13 : FlushCache(true);
333 13 : Mclose(d_map);
334 26 : }
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 9 : CPLErr PCRasterDataset::GetGeoTransform(GDALGeoTransform >) const
346 : {
347 : // x = west + nrCols * cellsize
348 9 : gt[0] = d_west;
349 9 : gt[1] = d_cellSize;
350 9 : gt[2] = 0.0;
351 :
352 : // y = north + nrRows * -cellsize
353 9 : gt[3] = d_north;
354 9 : gt[4] = 0.0;
355 9 : gt[5] = -1.0 * d_cellSize;
356 :
357 9 : return CE_None;
358 : }
359 :
360 : //! Returns the map handle.
361 : /*!
362 : \return Map handle.
363 : */
364 480 : MAP *PCRasterDataset::map() const
365 : {
366 480 : return d_map;
367 : }
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 433 : CSF_CR PCRasterDataset::cellRepresentation() const
376 : {
377 433 : return d_cellRepresentation;
378 : }
379 :
380 : //! Returns the value scale of the data.
381 : /*!
382 : \return Value scale
383 : \sa cellRepresentation()
384 : */
385 20 : CSF_VS PCRasterDataset::valueScale() const
386 : {
387 20 : return d_valueScale;
388 : }
389 :
390 : //! Returns the value of the missing value.
391 : /*!
392 : \return Missing value
393 : */
394 457 : double PCRasterDataset::defaultNoDataValue() const
395 : {
396 457 : return d_defaultNoDataValue;
397 : }
398 :
399 33 : GDALDataset *PCRasterDataset::create(const char *filename, int nr_cols,
400 : int nr_rows, int nrBands,
401 : GDALDataType gdalType,
402 : CSLConstList papszParamList)
403 : {
404 : // Checks
405 33 : if (nrBands != 1)
406 : {
407 18 : CPLError(CE_Failure, CPLE_NotSupported,
408 : "PCRaster driver : "
409 : "attempt to create dataset with too many bands (%d); "
410 : "must be 1 band.\n",
411 : nrBands);
412 18 : return nullptr;
413 : }
414 :
415 15 : const int row_col_max = INT4_MAX - 1;
416 15 : if (nr_cols > row_col_max)
417 : {
418 0 : CPLError(CE_Failure, CPLE_NotSupported,
419 : "PCRaster driver : "
420 : "attempt to create dataset with too many columns (%d); "
421 : "must be smaller than %d.",
422 : nr_cols, row_col_max);
423 0 : return nullptr;
424 : }
425 :
426 15 : if (nr_rows > row_col_max)
427 : {
428 0 : CPLError(CE_Failure, CPLE_NotSupported,
429 : "PCRaster driver : "
430 : "attempt to create dataset with too many rows (%d); "
431 : "must be smaller than %d.",
432 : nr_rows, row_col_max);
433 0 : return nullptr;
434 : }
435 :
436 15 : if (gdalType != GDT_UInt8 && gdalType != GDT_Int32 &&
437 : gdalType != GDT_Float32)
438 : {
439 11 : CPLError(CE_Failure, CPLE_AppDefined,
440 : "PCRaster driver: "
441 : "attempt to create dataset with an illegal data type (%s); "
442 : "use either Byte, Int32 or Float32.",
443 : GDALGetDataTypeName(gdalType));
444 11 : return nullptr;
445 : }
446 :
447 : // value scale must be specified by the user,
448 : // determines cell representation
449 : const char *valueScale =
450 4 : CSLFetchNameValue(papszParamList, "PCRASTER_VALUESCALE");
451 :
452 4 : if (valueScale == nullptr)
453 : {
454 3 : CPLError(CE_Failure, CPLE_AppDefined,
455 : "PCRaster driver: value scale can not be determined; "
456 : "specify PCRASTER_VALUESCALE.");
457 3 : return nullptr;
458 : }
459 :
460 1 : CSF_VS csf_value_scale = string2ValueScale(valueScale);
461 :
462 1 : if (csf_value_scale == VS_UNDEFINED)
463 : {
464 0 : CPLError(CE_Failure, CPLE_AppDefined,
465 : "PCRaster driver: value scale can not be determined (%s); "
466 : "use either VS_BOOLEAN, VS_NOMINAL, VS_ORDINAL, VS_SCALAR, "
467 : "VS_DIRECTION, VS_LDD",
468 : valueScale);
469 0 : return nullptr;
470 : }
471 :
472 : CSF_CR csf_cell_representation =
473 1 : GDALType2CellRepresentation(gdalType, false);
474 :
475 : // default values
476 1 : REAL8 west = 0.0;
477 1 : REAL8 north = 0.0;
478 1 : REAL8 length = 1.0;
479 1 : REAL8 angle = 0.0;
480 1 : CSF_PT projection = PT_YDECT2B;
481 :
482 : // Create a new raster
483 1 : MAP *map = Rcreate(filename, nr_rows, nr_cols, csf_cell_representation,
484 : csf_value_scale, projection, west, north, angle, length);
485 :
486 1 : if (!map)
487 : {
488 0 : CPLError(CE_Failure, CPLE_OpenFailed,
489 : "PCRaster driver: Unable to create raster %s", filename);
490 0 : return nullptr;
491 : }
492 :
493 1 : Mclose(map);
494 1 : map = nullptr;
495 :
496 1 : return GDALDataset::Open(filename, GDAL_OF_RASTER | GDAL_OF_UPDATE);
497 : }
498 :
499 0 : CPLErr PCRasterDataset::SetGeoTransform(const GDALGeoTransform >)
500 : {
501 0 : if ((gt[2] != 0.0) || (gt[4] != 0.0))
502 : {
503 0 : CPLError(CE_Failure, CPLE_NotSupported,
504 : "PCRaster driver: "
505 : "rotated geogtations are not supported.");
506 0 : return CE_Failure;
507 : }
508 :
509 0 : if (gt[1] != gt[5] * -1.0)
510 : {
511 0 : CPLError(CE_Failure, CPLE_NotSupported,
512 : "PCRaster driver: "
513 : "only the same width and height for cells is supported.");
514 0 : return CE_Failure;
515 : }
516 :
517 0 : d_west = gt[0];
518 0 : d_north = gt[3];
519 0 : d_cellSize = gt[1];
520 0 : d_location_changed = true;
521 :
522 0 : return CE_None;
523 : }
524 :
525 20 : bool PCRasterDataset::location_changed() const
526 : {
527 20 : return d_location_changed;
528 : }
|
