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