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