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