1 : /******************************************************************************
2 : * $Id: pcrasterdataset.cpp 16442 2009-03-01 16:46:45Z 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 16442 2009-03-01 16:46:45Z 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 9544 : GDALDataset* PCRasterDataset::open(
79 : GDALOpenInfo* info)
80 : {
81 9544 : PCRasterDataset* dataset = 0;
82 :
83 9544 : 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 9544 : 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 9544 : 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 16 : 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 16 : int nrBands = source->GetRasterCount();
135 16 : 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 11 : GDALRasterBand* raster = source->GetRasterBand(1);
142 :
143 : // Create PCRaster raster. Determine properties of raster to create.
144 11 : size_t nrRows = raster->GetYSize();
145 11 : size_t nrCols = raster->GetXSize();
146 11 : std::string string;
147 :
148 : // The in-file type of the cells.
149 : CSF_CR fileCellRepresentation = GDALType2CellRepresentation(
150 11 : raster->GetRasterDataType(), false);
151 :
152 11 : 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 7 : CSF_VS valueScale = VS_UNDEFINED;
160 7 : if(source->GetMetadataItem("PCRASTER_VALUESCALE")) {
161 0 : string = source->GetMetadataItem("PCRASTER_VALUESCALE");
162 : }
163 :
164 : valueScale = !string.empty()
165 : ? string2ValueScale(string)
166 7 : : GDALType2ValueScale(raster->GetRasterDataType());
167 :
168 7 : 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 7 : CSF_PT const projection = PT_YDECT2B;
175 7 : REAL8 const angle = 0.0;
176 7 : REAL8 west = 0.0;
177 7 : REAL8 north = 0.0;
178 7 : REAL8 cellSize = 1.0;
179 :
180 : double transform[6];
181 7 : if(source->GetGeoTransform(transform) == CE_None) {
182 7 : if(transform[2] == 0.0 && transform[4] == 0.0) {
183 7 : west = static_cast<REAL8>(transform[0]);
184 7 : north = static_cast<REAL8>(transform[3]);
185 7 : cellSize = static_cast<REAL8>(transform[1]);
186 : }
187 : }
188 :
189 : // The in-memory type of the cells.
190 7 : CSF_CR appCellRepresentation = CR_UNDEFINED;
191 : appCellRepresentation = GDALType2CellRepresentation(
192 7 : raster->GetRasterDataType(), true);
193 :
194 7 : 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 7 : valueScale = fitValueScale(valueScale, appCellRepresentation);
203 :
204 : // Create a raster with the in file cell representation.
205 : MAP* map = Rcreate(filename, nrRows, nrCols, fileCellRepresentation,
206 7 : valueScale, projection, west, north, angle, cellSize);
207 :
208 7 : if(!map) {
209 : CPLError(CE_Failure, CPLE_OpenFailed,
210 0 : "PCRaster driver: Unable to create raster %s", filename);
211 0 : 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 0 : free(buffer);
249 : CPLError(CE_Failure, CPLE_FileIO,
250 0 : "PCRaster driver: Error reading from source raster");
251 : }
252 :
253 : // Upon reading values are converted to the
254 : // right data type. This includes the missing value. If the source
255 : // value cannot be represented in the target data type it is set to a
256 : // missing value.
257 :
258 40 : if(hasMissingValue) {
259 0 : alterToStdMV(buffer, nrCols, appCellRepresentation, missingValue);
260 : }
261 :
262 40 : if(valueScale == VS_BOOLEAN) {
263 10 : castValuesToBooleanRange(buffer, nrCols, appCellRepresentation);
264 : }
265 :
266 : // Write row in target.
267 40 : RputRow(map, row, buffer);
268 :
269 40 : if(!progress((row + 1) / (static_cast<double>(nrRows)), 0, progressData)) {
270 0 : free(buffer);
271 : CPLError(CE_Failure, CPLE_UserInterrupt,
272 0 : "PCRaster driver: User terminated CreateCopy()");
273 : }
274 : }
275 :
276 4 : Mclose(map);
277 4 : map = 0;
278 :
279 4 : free(buffer);
280 4 : buffer = 0;
281 :
282 4 : if( errorCode != CE_None )
283 0 : return NULL;
284 :
285 : /* -------------------------------------------------------------------- */
286 : /* Re-open dataset, and copy any auxilary pam information. */
287 : /* -------------------------------------------------------------------- */
288 : GDALPamDataset *poDS = (GDALPamDataset *)
289 4 : GDALOpen( filename, GA_Update );
290 :
291 4 : if( poDS )
292 4 : poDS->CloneInfo( source, GCIF_PAM_DEFAULT );
293 :
294 4 : return poDS;
295 : }
296 :
297 :
298 :
299 : //------------------------------------------------------------------------------
300 : // DEFINITION OF PCRDATASET MEMBERS
301 : //------------------------------------------------------------------------------
302 :
303 : //! Constructor.
304 : /*!
305 : \param map PCRaster map handle. It is ours to close.
306 : */
307 6 : PCRasterDataset::PCRasterDataset(
308 : MAP* map)
309 :
310 : : GDALPamDataset(),
311 6 : d_map(map), d_west(0.0), d_north(0.0), d_cellSize(0.0)
312 :
313 : {
314 : // Read header info.
315 6 : nRasterXSize = RgetNrCols(d_map);
316 6 : nRasterYSize = RgetNrRows(d_map);
317 6 : d_west = static_cast<double>(RgetXUL(d_map));
318 6 : d_north = static_cast<double>(RgetYUL(d_map));
319 6 : d_cellSize = static_cast<double>(RgetCellSize(d_map));
320 6 : d_cellRepresentation = RgetUseCellRepr(d_map);
321 : CPLAssert(d_cellRepresentation != CR_UNDEFINED);
322 6 : d_valueScale = RgetValueScale(d_map);
323 : CPLAssert(d_valueScale != VS_UNDEFINED);
324 6 : d_missingValue = ::missingValue(d_cellRepresentation);
325 :
326 : // Create band information objects.
327 6 : nBands = 1;
328 6 : SetBand(1, new PCRasterRasterBand(this));
329 :
330 : SetMetadataItem("PCRASTER_VALUESCALE", valueScale2String(
331 6 : d_valueScale).c_str());
332 6 : }
333 :
334 :
335 :
336 : //! Destructor.
337 : /*!
338 : \warning The map given in the constructor is closed.
339 : */
340 12 : PCRasterDataset::~PCRasterDataset()
341 : {
342 6 : FlushCache();
343 6 : Mclose(d_map);
344 12 : }
345 :
346 :
347 :
348 : //! Sets projections info.
349 : /*!
350 : \param transform Array to fill.
351 :
352 : CSF 2.0 supports the notion of y coordinates which increase from north to
353 : south. Support for this has been dropped and applications reading PCRaster
354 : rasters will treat or already treat y coordinates as increasing from south
355 : to north only.
356 : */
357 5 : CPLErr PCRasterDataset::GetGeoTransform(double* transform)
358 : {
359 : // x = west + nrCols * cellsize
360 5 : transform[0] = d_west;
361 5 : transform[1] = d_cellSize;
362 5 : transform[2] = 0.0;
363 :
364 : // y = north + nrRows * -cellsize
365 5 : transform[3] = d_north;
366 5 : transform[4] = 0.0;
367 5 : transform[5] = -1.0 * d_cellSize;
368 :
369 5 : return CE_None;
370 : }
371 :
372 :
373 :
374 : //! Returns the map handle.
375 : /*!
376 : \return Map handle.
377 : */
378 100 : MAP* PCRasterDataset::map() const
379 : {
380 100 : return d_map;
381 : }
382 :
383 :
384 :
385 : //! Returns the in-app cell representation.
386 : /*!
387 : \return cell representation
388 : \warning This might not be the same representation as use to store the values in the file.
389 : \sa valueScale()
390 : */
391 106 : CSF_CR PCRasterDataset::cellRepresentation() const
392 : {
393 106 : return d_cellRepresentation;
394 : }
395 :
396 :
397 :
398 : //! Returns the value scale of the data.
399 : /*!
400 : \return Value scale
401 : \sa cellRepresentation()
402 : */
403 0 : CSF_VS PCRasterDataset::valueScale() const
404 : {
405 0 : return d_valueScale;
406 : }
407 :
408 :
409 :
410 : //! Returns the value of the missing value.
411 : /*!
412 : \return Missing value
413 : */
414 101 : double PCRasterDataset::missingValue() const
415 : {
416 101 : return d_missingValue;
417 1140 : }
418 :
419 :
420 :
421 : //------------------------------------------------------------------------------
422 : // DEFINITION OF FREE OPERATORS
423 : //------------------------------------------------------------------------------
424 :
425 :
426 :
427 : //------------------------------------------------------------------------------
428 : // DEFINITION OF FREE FUNCTIONS
429 : //------------------------------------------------------------------------------
430 :
431 :
432 :
|
