1 : /******************************************************************************
2 : * $Id: pcrasterdataset.cpp 22609 2011-06-28 21:01:48Z rouault $
3 : *
4 : * Project: PCRaster Integration
5 : * Purpose: PCRaster CSF 2.0 raster file driver
6 : * Author: Kor de Jong, k.dejong at geog.uu.nl
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Kor de Jong
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 : #include "gdal_pam.h"
31 : #include "cpl_string.h"
32 :
33 : CPL_CVSID("$Id: pcrasterdataset.cpp 22609 2011-06-28 21:01:48Z rouault $");
34 :
35 : #ifndef INCLUDED_PCRASTERDATASET
36 : #include "pcrasterdataset.h"
37 : #define INCLUDED_PCRASTERDATASET
38 : #endif
39 :
40 : #ifndef INCLUDED_IOSTREAM
41 : #include <iostream>
42 : #define INCLUDED_IOSTREAM
43 : #endif
44 :
45 : // PCRaster library headers.
46 :
47 : // Module headers.
48 : #ifndef INCLUDED_PCRASTERRASTERBAND
49 : #include "pcrasterrasterband.h"
50 : #define INCLUDED_PCRASTERRASTERBAND
51 : #endif
52 :
53 : #ifndef INCLUDED_PCRASTERUTIL
54 : #include "pcrasterutil.h"
55 : #define INCLUDED_PCRASTERUTIL
56 : #endif
57 :
58 :
59 :
60 : /*!
61 : \file
62 : This file contains the implementation of the PCRasterDataset class.
63 : */
64 :
65 :
66 :
67 : //------------------------------------------------------------------------------
68 : // DEFINITION OF STATIC PCRDATASET MEMBERS
69 : //------------------------------------------------------------------------------
70 :
71 : //! Tries to open the file described by \a info.
72 : /*!
73 : \param info Object with information about the dataset to open.
74 : \return Pointer to newly allocated GDALDataset or 0.
75 :
76 : Returns 0 if the file could not be opened.
77 : */
78 13562 : GDALDataset* PCRasterDataset::open(
79 : GDALOpenInfo* info)
80 : {
81 13562 : PCRasterDataset* dataset = 0;
82 :
83 13562 : if(info->fp && info->nHeaderBytes >= static_cast<int>(CSF_SIZE_SIG) &&
84 : strncmp((char*)info->pabyHeader, CSF_SIG, CSF_SIZE_SIG) == 0) {
85 : MOPEN_PERM mode = info->eAccess == GA_Update
86 : ? M_READ_WRITE
87 6 : : M_READ;
88 :
89 6 : MAP* map = mapOpen(info->pszFilename, mode);
90 :
91 6 : if(map) {
92 6 : dataset = new PCRasterDataset(map);
93 : }
94 : }
95 :
96 : /* -------------------------------------------------------------------- */
97 : /* Initialize any PAM information and overviews. */
98 : /* -------------------------------------------------------------------- */
99 13562 : if( dataset )
100 : {
101 6 : dataset->SetDescription( info->pszFilename );
102 6 : dataset->TryLoadXML();
103 :
104 6 : dataset->oOvManager.Initialize( dataset, info->pszFilename );
105 : }
106 :
107 13562 : return dataset;
108 : }
109 :
110 :
111 :
112 : //! Writes a raster to \a filename as a PCRaster raster file.
113 : /*!
114 : \warning The source raster must have only 1 band. Currently, the values in
115 : the source raster must be stored in one of the supported cell
116 : representations (CR_UINT1, CR_INT4, CR_REAL4, CR_REAL8).
117 :
118 : The meta data item PCRASTER_VALUESCALE will be checked to see what value
119 : scale to use. Otherwise a value scale is determined using
120 : GDALType2ValueScale(GDALDataType).
121 :
122 : This function always writes rasters using CR_UINT1, CR_INT4 or CR_REAL4
123 : cell representations.
124 : */
125 22 : GDALDataset* PCRasterDataset::createCopy(
126 : char const* filename,
127 : GDALDataset* source,
128 : int strict,
129 : char** options,
130 : GDALProgressFunc progress,
131 : void* progressData)
132 : {
133 : // Checks.
134 22 : int nrBands = source->GetRasterCount();
135 22 : if(nrBands != 1) {
136 : CPLError(CE_Failure, CPLE_NotSupported,
137 5 : "PCRaster driver: Too many bands ('%d'): must be 1 band", nrBands);
138 5 : return 0;
139 : }
140 :
141 17 : GDALRasterBand* raster = source->GetRasterBand(1);
142 :
143 : // Create PCRaster raster. Determine properties of raster to create.
144 17 : size_t nrRows = raster->GetYSize();
145 17 : size_t nrCols = raster->GetXSize();
146 17 : std::string string;
147 :
148 : // The in-file type of the cells.
149 : CSF_CR fileCellRepresentation = GDALType2CellRepresentation(
150 17 : raster->GetRasterDataType(), false);
151 :
152 17 : if(fileCellRepresentation == CR_UNDEFINED) {
153 : CPLError(CE_Failure, CPLE_NotSupported,
154 4 : "PCRaster driver: Cannot determine a valid cell representation");
155 4 : return 0;
156 : }
157 :
158 : // The value scale of the values.
159 13 : CSF_VS valueScale = VS_UNDEFINED;
160 13 : if(source->GetMetadataItem("PCRASTER_VALUESCALE")) {
161 0 : string = source->GetMetadataItem("PCRASTER_VALUESCALE");
162 : }
163 :
164 : valueScale = !string.empty()
165 : ? string2ValueScale(string)
166 13 : : GDALType2ValueScale(raster->GetRasterDataType());
167 :
168 13 : if(valueScale == VS_UNDEFINED) {
169 : CPLError(CE_Failure, CPLE_NotSupported,
170 0 : "PCRaster driver: Cannot determine a valid value scale");
171 0 : return 0;
172 : }
173 :
174 13 : CSF_PT const projection = PT_YDECT2B;
175 13 : REAL8 const angle = 0.0;
176 13 : REAL8 west = 0.0;
177 13 : REAL8 north = 0.0;
178 13 : REAL8 cellSize = 1.0;
179 :
180 : double transform[6];
181 13 : if(source->GetGeoTransform(transform) == CE_None) {
182 13 : if(transform[2] == 0.0 && transform[4] == 0.0) {
183 13 : west = static_cast<REAL8>(transform[0]);
184 13 : north = static_cast<REAL8>(transform[3]);
185 13 : cellSize = static_cast<REAL8>(transform[1]);
186 : }
187 : }
188 :
189 : // The in-memory type of the cells.
190 13 : CSF_CR appCellRepresentation = CR_UNDEFINED;
191 : appCellRepresentation = GDALType2CellRepresentation(
192 13 : raster->GetRasterDataType(), true);
193 :
194 13 : if(appCellRepresentation == CR_UNDEFINED) {
195 : CPLError(CE_Failure, CPLE_NotSupported,
196 0 : "PCRaster driver: Cannot determine a valid cell representation");
197 0 : return 0;
198 : }
199 :
200 : // Check whether value scale fits the cell representation. Adjust when
201 : // needed.
202 13 : valueScale = fitValueScale(valueScale, appCellRepresentation);
203 :
204 : // Create a raster with the in file cell representation.
205 : MAP* map = Rcreate(filename, nrRows, nrCols, fileCellRepresentation,
206 13 : valueScale, projection, west, north, angle, cellSize);
207 :
208 13 : if(!map) {
209 : CPLError(CE_Failure, CPLE_OpenFailed,
210 6 : "PCRaster driver: Unable to create raster %s", filename);
211 6 : return 0;
212 : }
213 :
214 : // Try to convert in app cell representation to the cell representation
215 : // of the file.
216 7 : if(RuseAs(map, appCellRepresentation)) {
217 : CPLError(CE_Failure, CPLE_NotSupported,
218 3 : "PCRaster driver: Cannot convert cells: %s", MstrError());
219 3 : Mclose(map);
220 3 : return 0;
221 : }
222 :
223 : int hasMissingValue;
224 4 : double missingValue = raster->GetNoDataValue(&hasMissingValue);
225 :
226 : // This is needed to get my (KDJ) unit tests running.
227 : // I am still uncertain why this is needed. If the input raster has float32
228 : // values and the output int32, than the missing value in the dataset object
229 : // is not updated like the values are.
230 4 : if(missingValue == ::missingValue(CR_REAL4) &&
231 : fileCellRepresentation == CR_INT4) {
232 0 : missingValue = ::missingValue(fileCellRepresentation);
233 : }
234 :
235 : // TODO conversie van INT2 naar INT4 ondersteunen. zie ruseas.c regel 503.
236 : // conversie op r 159.
237 :
238 : // Create buffer for one row of values.
239 4 : void* buffer = Rmalloc(map, nrCols);
240 :
241 : // Copy values from source to target.
242 4 : CPLErr errorCode = CE_None;
243 44 : for(size_t row = 0; errorCode == CE_None && row < nrRows; ++row) {
244 :
245 : // Get row from source.
246 40 : if(raster->RasterIO(GF_Read, 0, row, nrCols, 1, buffer, nrCols, 1,
247 : raster->GetRasterDataType(), 0, 0) != CE_None) {
248 : CPLError(CE_Failure, CPLE_FileIO,
249 0 : "PCRaster driver: Error reading from source raster");
250 0 : errorCode = CE_Failure;
251 0 : break;
252 : }
253 :
254 : // Upon reading values are converted to the
255 : // right data type. This includes the missing value. If the source
256 : // value cannot be represented in the target data type it is set to a
257 : // missing value.
258 :
259 40 : if(hasMissingValue) {
260 0 : alterToStdMV(buffer, nrCols, appCellRepresentation, missingValue);
261 : }
262 :
263 40 : if(valueScale == VS_BOOLEAN) {
264 10 : castValuesToBooleanRange(buffer, nrCols, appCellRepresentation);
265 : }
266 :
267 : // Write row in target.
268 40 : RputRow(map, row, buffer);
269 :
270 40 : if(!progress((row + 1) / (static_cast<double>(nrRows)), 0, progressData)) {
271 : CPLError(CE_Failure, CPLE_UserInterrupt,
272 0 : "PCRaster driver: User terminated CreateCopy()");
273 0 : errorCode = CE_Failure;
274 0 : break;
275 : }
276 : }
277 :
278 4 : Mclose(map);
279 4 : map = 0;
280 :
281 4 : free(buffer);
282 4 : buffer = 0;
283 :
284 4 : if( errorCode != CE_None )
285 0 : return NULL;
286 :
287 : /* -------------------------------------------------------------------- */
288 : /* Re-open dataset, and copy any auxilary pam information. */
289 : /* -------------------------------------------------------------------- */
290 : GDALPamDataset *poDS = (GDALPamDataset *)
291 4 : GDALOpen( filename, GA_Update );
292 :
293 4 : if( poDS )
294 4 : poDS->CloneInfo( source, GCIF_PAM_DEFAULT );
295 :
296 4 : return poDS;
297 : }
298 :
299 :
300 :
301 : //------------------------------------------------------------------------------
302 : // DEFINITION OF PCRDATASET MEMBERS
303 : //------------------------------------------------------------------------------
304 :
305 : //! Constructor.
306 : /*!
307 : \param map PCRaster map handle. It is ours to close.
308 : */
309 6 : PCRasterDataset::PCRasterDataset(
310 : MAP* map)
311 :
312 : : GDALPamDataset(),
313 6 : d_map(map), d_west(0.0), d_north(0.0), d_cellSize(0.0)
314 :
315 : {
316 : // Read header info.
317 6 : nRasterXSize = RgetNrCols(d_map);
318 6 : nRasterYSize = RgetNrRows(d_map);
319 6 : d_west = static_cast<double>(RgetXUL(d_map));
320 6 : d_north = static_cast<double>(RgetYUL(d_map));
321 6 : d_cellSize = static_cast<double>(RgetCellSize(d_map));
322 6 : d_cellRepresentation = RgetUseCellRepr(d_map);
323 6 : CPLAssert(d_cellRepresentation != CR_UNDEFINED);
324 6 : d_valueScale = RgetValueScale(d_map);
325 6 : CPLAssert(d_valueScale != VS_UNDEFINED);
326 6 : d_missingValue = ::missingValue(d_cellRepresentation);
327 :
328 : // Create band information objects.
329 6 : nBands = 1;
330 6 : SetBand(1, new PCRasterRasterBand(this));
331 :
332 : SetMetadataItem("PCRASTER_VALUESCALE", valueScale2String(
333 6 : d_valueScale).c_str());
334 6 : }
335 :
336 :
337 :
338 : //! Destructor.
339 : /*!
340 : \warning The map given in the constructor is closed.
341 : */
342 6 : PCRasterDataset::~PCRasterDataset()
343 : {
344 6 : FlushCache();
345 6 : Mclose(d_map);
346 6 : }
347 :
348 :
349 :
350 : //! Sets projections info.
351 : /*!
352 : \param transform Array to fill.
353 :
354 : CSF 2.0 supports the notion of y coordinates which increase from north to
355 : south. Support for this has been dropped and applications reading PCRaster
356 : rasters will treat or already treat y coordinates as increasing from south
357 : to north only.
358 : */
359 5 : CPLErr PCRasterDataset::GetGeoTransform(double* transform)
360 : {
361 : // x = west + nrCols * cellsize
362 5 : transform[0] = d_west;
363 5 : transform[1] = d_cellSize;
364 5 : transform[2] = 0.0;
365 :
366 : // y = north + nrRows * -cellsize
367 5 : transform[3] = d_north;
368 5 : transform[4] = 0.0;
369 5 : transform[5] = -1.0 * d_cellSize;
370 :
371 5 : return CE_None;
372 : }
373 :
374 :
375 :
376 : //! Returns the map handle.
377 : /*!
378 : \return Map handle.
379 : */
380 100 : MAP* PCRasterDataset::map() const
381 : {
382 100 : return d_map;
383 : }
384 :
385 :
386 :
387 : //! Returns the in-app cell representation.
388 : /*!
389 : \return cell representation
390 : \warning This might not be the same representation as use to store the values in the file.
391 : \sa valueScale()
392 : */
393 106 : CSF_CR PCRasterDataset::cellRepresentation() const
394 : {
395 106 : return d_cellRepresentation;
396 : }
397 :
398 :
399 :
400 : //! Returns the value scale of the data.
401 : /*!
402 : \return Value scale
403 : \sa cellRepresentation()
404 : */
405 0 : CSF_VS PCRasterDataset::valueScale() const
406 : {
407 0 : return d_valueScale;
408 : }
409 :
410 :
411 :
412 : //! Returns the value of the missing value.
413 : /*!
414 : \return Missing value
415 : */
416 101 : double PCRasterDataset::missingValue() const
417 : {
418 101 : return d_missingValue;
419 2139 : }
420 :
421 :
422 :
423 : //------------------------------------------------------------------------------
424 : // DEFINITION OF FREE OPERATORS
425 : //------------------------------------------------------------------------------
426 :
427 :
428 :
429 : //------------------------------------------------------------------------------
430 : // DEFINITION OF FREE FUNCTIONS
431 : //------------------------------------------------------------------------------
432 :
433 :
434 :
|
