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