1 : /******************************************************************************
2 : * $Id: rdataset.cpp 17718 2009-09-30 20:15:36Z warmerdam $
3 : *
4 : * Project: R Format Driver
5 : * Purpose: Read/write R stats package object format.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2009, Frank Warmerdam <warmerdam@pobox.com>
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 : #include "../raw/rawdataset.h"
33 :
34 : CPL_CVSID("$Id: rdataset.cpp 17718 2009-09-30 20:15:36Z warmerdam $");
35 :
36 : CPL_C_START
37 : void GDALRegister_R(void);
38 : CPL_C_END
39 :
40 : GDALDataset *
41 : RCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
42 : int bStrict, char ** papszOptions,
43 : GDALProgressFunc pfnProgress, void * pProgressData );
44 :
45 : #define R_NILSXP 0
46 : #define R_LISTSXP 2
47 : #define R_CHARSXP 9
48 : #define R_INTSXP 13
49 : #define R_REALSXP 14
50 : #define R_STRSXP 16
51 :
52 : /************************************************************************/
53 : /* ==================================================================== */
54 : /* RDataset */
55 : /* ==================================================================== */
56 : /************************************************************************/
57 :
58 : class RDataset : public GDALPamDataset
59 : {
60 : friend class RRasterBand;
61 : FILE *fp;
62 : int bASCII;
63 : CPLString osLastStringRead;
64 :
65 : vsi_l_offset nStartOfData;
66 :
67 : double *padfMatrixValues;
68 :
69 : const char *ASCIIFGets();
70 : int ReadInteger();
71 : double ReadFloat();
72 : const char *ReadString();
73 : int ReadPair( CPLString &osItemName, int &nItemType );
74 :
75 : public:
76 : RDataset();
77 : ~RDataset();
78 :
79 : static GDALDataset *Open( GDALOpenInfo * );
80 : static int Identify( GDALOpenInfo * );
81 : };
82 :
83 : /************************************************************************/
84 : /* ==================================================================== */
85 : /* RRasterBand */
86 : /* ==================================================================== */
87 : /************************************************************************/
88 :
89 : class RRasterBand : public GDALPamRasterBand
90 : {
91 : friend class RDataset;
92 :
93 : const double *padfMatrixValues;
94 :
95 : public:
96 :
97 : RRasterBand( RDataset *, int, const double * );
98 : ~RRasterBand();
99 :
100 : virtual CPLErr IReadBlock( int, int, void * );
101 : };
102 :
103 : /************************************************************************/
104 : /* RRasterBand() */
105 : /************************************************************************/
106 :
107 5 : RRasterBand::RRasterBand( RDataset *poDS, int nBand,
108 5 : const double *padfMatrixValues )
109 : {
110 5 : this->poDS = poDS;
111 5 : this->nBand = nBand;
112 5 : this->padfMatrixValues = padfMatrixValues;
113 :
114 5 : eDataType = GDT_Float64;
115 :
116 5 : nBlockXSize = poDS->nRasterXSize;
117 5 : nBlockYSize = 1;
118 5 : }
119 :
120 : /************************************************************************/
121 : /* ~RRasterBand() */
122 : /************************************************************************/
123 :
124 10 : RRasterBand::~RRasterBand()
125 : {
126 10 : }
127 :
128 : /************************************************************************/
129 : /* IReadBlock() */
130 : /************************************************************************/
131 :
132 45 : CPLErr RRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
133 : void * pImage )
134 : {
135 : memcpy( pImage, padfMatrixValues + nBlockYOff * nBlockXSize,
136 45 : nBlockXSize * 8 );
137 45 : return CE_None;
138 : }
139 :
140 : /************************************************************************/
141 : /* ==================================================================== */
142 : /* RDataset() */
143 : /* ==================================================================== */
144 : /************************************************************************/
145 :
146 : /************************************************************************/
147 : /* RDataset() */
148 : /************************************************************************/
149 :
150 8 : RDataset::RDataset()
151 : {
152 8 : fp = NULL;
153 8 : padfMatrixValues = NULL;
154 8 : }
155 :
156 : /************************************************************************/
157 : /* ~RDataset() */
158 : /************************************************************************/
159 :
160 16 : RDataset::~RDataset()
161 : {
162 8 : FlushCache();
163 8 : CPLFree(padfMatrixValues);
164 :
165 8 : if( fp )
166 8 : VSIFCloseL( fp );
167 16 : }
168 :
169 : /************************************************************************/
170 : /* ASCIIFGets() */
171 : /* */
172 : /* Fetch one line from an ASCII source into osLastStringRead. */
173 : /************************************************************************/
174 :
175 1324 : const char *RDataset::ASCIIFGets()
176 :
177 : {
178 : char chNextChar;
179 :
180 1324 : osLastStringRead.resize(0);
181 :
182 5109 : do
183 : {
184 5109 : chNextChar = '\n';
185 5109 : VSIFReadL( &chNextChar, 1, 1, fp );
186 5109 : if( chNextChar != '\n' )
187 3785 : osLastStringRead += chNextChar;
188 : } while( chNextChar != '\n' && chNextChar != '\0' );
189 :
190 1324 : return osLastStringRead;
191 : }
192 :
193 : /************************************************************************/
194 : /* ReadInteger() */
195 : /************************************************************************/
196 :
197 152 : int RDataset::ReadInteger()
198 :
199 : {
200 152 : if( bASCII )
201 : {
202 76 : return atoi(ASCIIFGets());
203 : }
204 : else
205 : {
206 : GInt32 nValue;
207 :
208 76 : if( VSIFReadL( &nValue, 4, 1, fp ) != 1 )
209 0 : return -1;
210 76 : CPL_MSBPTR32( &nValue );
211 :
212 76 : return nValue;
213 : }
214 : }
215 :
216 : /************************************************************************/
217 : /* ReadFloat() */
218 : /************************************************************************/
219 :
220 1240 : double RDataset::ReadFloat()
221 :
222 : {
223 1240 : if( bASCII )
224 : {
225 1240 : return atof(ASCIIFGets());
226 : }
227 : else
228 : {
229 : double dfValue;
230 :
231 0 : if( VSIFReadL( &dfValue, 8, 1, fp ) != 1 )
232 0 : return -1;
233 0 : CPL_MSBPTR64( &dfValue );
234 :
235 0 : return dfValue;
236 : }
237 : }
238 :
239 : /************************************************************************/
240 : /* ReadString() */
241 : /************************************************************************/
242 :
243 16 : const char *RDataset::ReadString()
244 :
245 : {
246 16 : if( ReadInteger() % 256 != R_CHARSXP )
247 : {
248 0 : osLastStringRead = "";
249 0 : return "";
250 : }
251 :
252 16 : size_t nLen = ReadInteger();
253 :
254 16 : char *pachWrkBuf = (char *) VSIMalloc(nLen);
255 16 : if (pachWrkBuf == NULL)
256 : {
257 0 : osLastStringRead = "";
258 0 : return "";
259 : }
260 16 : if( VSIFReadL( pachWrkBuf, 1, nLen, fp ) != nLen )
261 : {
262 0 : osLastStringRead = "";
263 0 : CPLFree( pachWrkBuf );
264 0 : return "";
265 : }
266 :
267 16 : if( bASCII )
268 : {
269 : /* suck up newline and any extra junk */
270 8 : ASCIIFGets();
271 : }
272 :
273 16 : osLastStringRead.assign( pachWrkBuf, nLen );
274 16 : CPLFree( pachWrkBuf );
275 :
276 16 : return osLastStringRead;
277 : }
278 :
279 : /************************************************************************/
280 : /* ReadPair() */
281 : /************************************************************************/
282 :
283 24 : int RDataset::ReadPair( CPLString &osObjName, int &nObjCode )
284 :
285 : {
286 24 : nObjCode = ReadInteger();
287 24 : if( nObjCode == 254 )
288 8 : return TRUE;
289 :
290 16 : if( (nObjCode % 256) != R_LISTSXP )
291 : {
292 : CPLError( CE_Failure, CPLE_OpenFailed,
293 0 : "Did not find expected object pair object." );
294 0 : return FALSE;
295 : }
296 :
297 16 : int nPairCount = ReadInteger();
298 16 : if( nPairCount != 1 )
299 : {
300 : CPLError( CE_Failure, CPLE_OpenFailed,
301 0 : "Did not find expected pair count of 1." );
302 0 : return FALSE;
303 : }
304 :
305 : /* -------------------------------------------------------------------- */
306 : /* Read the object name. */
307 : /* -------------------------------------------------------------------- */
308 16 : const char *pszName = ReadString();
309 16 : if( pszName == NULL || pszName[0] == '\0' )
310 0 : return FALSE;
311 :
312 16 : osObjName = pszName;
313 :
314 : /* -------------------------------------------------------------------- */
315 : /* Confirm that we have a numeric matrix object. */
316 : /* -------------------------------------------------------------------- */
317 16 : nObjCode = ReadInteger();
318 :
319 16 : return TRUE;
320 : }
321 :
322 : /************************************************************************/
323 : /* Identify() */
324 : /************************************************************************/
325 :
326 8784 : int RDataset::Identify( GDALOpenInfo *poOpenInfo )
327 : {
328 8784 : if( poOpenInfo->nHeaderBytes < 50 )
329 8425 : return FALSE;
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* If the extension is .rda and the file type is gzip */
333 : /* compressed we assume it is a gziped R binary file. */
334 : /* -------------------------------------------------------------------- */
335 359 : if( memcmp(poOpenInfo->pabyHeader,"\037\213\b",3) == 0
336 : && EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"rda") )
337 3 : return TRUE;
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Is this an ASCII or XDR binary R file? */
341 : /* -------------------------------------------------------------------- */
342 356 : if( !EQUALN((const char *)poOpenInfo->pabyHeader,"RDA2\nA\n",7)
343 : && !EQUALN((const char *)poOpenInfo->pabyHeader,"RDX2\nX\n",7) )
344 351 : return FALSE;
345 :
346 5 : return TRUE;
347 : }
348 :
349 : /************************************************************************/
350 : /* Open() */
351 : /************************************************************************/
352 :
353 1057 : GDALDataset *RDataset::Open( GDALOpenInfo * poOpenInfo )
354 : {
355 1057 : if( !Identify( poOpenInfo ) )
356 1049 : return NULL;
357 :
358 : /* -------------------------------------------------------------------- */
359 : /* Confirm the requested access is supported. */
360 : /* -------------------------------------------------------------------- */
361 8 : if( poOpenInfo->eAccess == GA_Update )
362 : {
363 : CPLError( CE_Failure, CPLE_NotSupported,
364 : "The R driver does not support update access to existing"
365 0 : " datasets.\n" );
366 0 : return NULL;
367 : }
368 :
369 : /* -------------------------------------------------------------------- */
370 : /* Do we need to route the file through the decompression */
371 : /* machinery? */
372 : /* -------------------------------------------------------------------- */
373 8 : CPLString osAdjustedFilename;
374 :
375 8 : if( memcmp(poOpenInfo->pabyHeader,"\037\213\b",3) == 0 )
376 3 : osAdjustedFilename = "/vsigzip/";
377 :
378 8 : osAdjustedFilename += poOpenInfo->pszFilename;
379 :
380 : /* -------------------------------------------------------------------- */
381 : /* Establish this as a dataset and open the file using VSI*L. */
382 : /* -------------------------------------------------------------------- */
383 8 : RDataset *poDS = new RDataset();
384 :
385 16 : poDS->fp = VSIFOpenL( osAdjustedFilename, "r" );
386 8 : if( poDS->fp == NULL )
387 : {
388 0 : delete poDS;
389 0 : return NULL;
390 : }
391 :
392 8 : poDS->bASCII = EQUALN((const char *)poOpenInfo->pabyHeader,"RDA2\nA\n",7);
393 :
394 : /* -------------------------------------------------------------------- */
395 : /* Confirm this is a version 2 file. */
396 : /* -------------------------------------------------------------------- */
397 8 : VSIFSeekL( poDS->fp, 7, SEEK_SET );
398 8 : if( poDS->ReadInteger() != R_LISTSXP )
399 : {
400 0 : delete poDS;
401 : CPLError( CE_Failure, CPLE_OpenFailed,
402 : "It appears %s is not a version 2 R object file after all!",
403 0 : poOpenInfo->pszFilename );
404 0 : return NULL;
405 : }
406 :
407 : /* -------------------------------------------------------------------- */
408 : /* Skip the version values. */
409 : /* -------------------------------------------------------------------- */
410 8 : poDS->ReadInteger();
411 8 : poDS->ReadInteger();
412 :
413 : /* -------------------------------------------------------------------- */
414 : /* Confirm we have a numeric vector object in a pairlist. */
415 : /* -------------------------------------------------------------------- */
416 8 : CPLString osObjName;
417 : int nObjCode;
418 :
419 8 : if( !poDS->ReadPair( osObjName, nObjCode ) )
420 : {
421 0 : delete poDS;
422 0 : return NULL;
423 : }
424 :
425 8 : if( nObjCode % 256 != R_REALSXP )
426 : {
427 0 : delete poDS;
428 : CPLError( CE_Failure, CPLE_OpenFailed,
429 0 : "Failed to find expected numeric vector object." );
430 0 : return NULL;
431 : }
432 :
433 8 : poDS->SetMetadataItem( "R_OBJECT_NAME", osObjName );
434 :
435 : /* -------------------------------------------------------------------- */
436 : /* Read the count. */
437 : /* -------------------------------------------------------------------- */
438 8 : int nValueCount = poDS->ReadInteger();
439 :
440 8 : poDS->nStartOfData = VSIFTellL( poDS->fp );
441 :
442 : /* -------------------------------------------------------------------- */
443 : /* Read/Skip ahead to attributes. */
444 : /* -------------------------------------------------------------------- */
445 8 : if( poDS->bASCII )
446 : {
447 4 : poDS->padfMatrixValues = (double*) VSIMalloc2( nValueCount, sizeof(double) );
448 4 : if (poDS->padfMatrixValues == NULL)
449 : {
450 : CPLError(CE_Failure, CPLE_AppDefined,
451 0 : "Cannot allocate %d doubles", nValueCount);
452 0 : delete poDS;
453 0 : return NULL;
454 : }
455 1244 : for( int iValue = 0; iValue < nValueCount; iValue++ )
456 1240 : poDS->padfMatrixValues[iValue] = poDS->ReadFloat();
457 : }
458 : else
459 : {
460 4 : VSIFSeekL( poDS->fp, 8 * nValueCount, SEEK_CUR );
461 : }
462 :
463 : /* -------------------------------------------------------------------- */
464 : /* Read pairs till we run out, trying to find a few items that */
465 : /* have special meaning to us. */
466 : /* -------------------------------------------------------------------- */
467 8 : poDS->nRasterXSize = poDS->nRasterYSize = 0;
468 8 : int nBandCount = 0;
469 :
470 24 : while( poDS->ReadPair( osObjName, nObjCode ) && nObjCode != 254 )
471 : {
472 8 : if( osObjName == "dim" && nObjCode % 256 == R_INTSXP )
473 : {
474 8 : int nCount = poDS->ReadInteger();
475 8 : if( nCount == 2 )
476 : {
477 0 : poDS->nRasterXSize = poDS->ReadInteger();
478 0 : poDS->nRasterYSize = poDS->ReadInteger();
479 0 : nBandCount = 1;
480 : }
481 8 : else if( nCount == 3 )
482 : {
483 8 : poDS->nRasterXSize = poDS->ReadInteger();
484 8 : poDS->nRasterYSize = poDS->ReadInteger();
485 8 : nBandCount = poDS->ReadInteger();
486 : }
487 : else
488 : {
489 : CPLError( CE_Failure, CPLE_AppDefined,
490 0 : "R 'dim' dimension wrong." );
491 0 : delete poDS;
492 0 : return NULL;
493 : }
494 : }
495 0 : else if( nObjCode % 256 == R_REALSXP )
496 : {
497 0 : int nCount = poDS->ReadInteger();
498 0 : while( nCount-- > 0 && !VSIFEofL(poDS->fp) )
499 0 : poDS->ReadFloat();
500 : }
501 0 : else if( nObjCode % 256 == R_INTSXP )
502 : {
503 0 : int nCount = poDS->ReadInteger();
504 0 : while( nCount-- > 0 && !VSIFEofL(poDS->fp) )
505 0 : poDS->ReadInteger();
506 : }
507 0 : else if( nObjCode % 256 == R_STRSXP )
508 : {
509 0 : int nCount = poDS->ReadInteger();
510 0 : while( nCount-- > 0 && !VSIFEofL(poDS->fp) )
511 0 : poDS->ReadString();
512 : }
513 0 : else if( nObjCode % 256 == R_CHARSXP )
514 : {
515 0 : poDS->ReadString();
516 : }
517 : }
518 :
519 8 : if( poDS->nRasterXSize == 0 )
520 : {
521 0 : delete poDS;
522 : CPLError( CE_Failure, CPLE_AppDefined,
523 0 : "Failed to find dim dimension information for R dataset." );
524 0 : return NULL;
525 : }
526 :
527 8 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
528 : !GDALCheckBandCount(nBandCount, TRUE))
529 : {
530 0 : delete poDS;
531 0 : return NULL;
532 : }
533 :
534 8 : if( nValueCount
535 : < ((GIntBig) nBandCount) * poDS->nRasterXSize * poDS->nRasterYSize )
536 : {
537 : CPLError( CE_Failure, CPLE_AppDefined,
538 0 : "Not enough pixel data." );
539 0 : delete poDS;
540 0 : return NULL;
541 : }
542 :
543 : /* -------------------------------------------------------------------- */
544 : /* Create the raster band object(s). */
545 : /* -------------------------------------------------------------------- */
546 18 : for( int iBand = 0; iBand < nBandCount; iBand++ )
547 : {
548 : GDALRasterBand *poBand;
549 :
550 10 : if( poDS->bASCII )
551 : poBand = new RRasterBand( poDS, iBand+1,
552 5 : poDS->padfMatrixValues + iBand * poDS->nRasterXSize * poDS->nRasterYSize );
553 : else
554 : poBand = new RawRasterBand( poDS, iBand+1, poDS->fp,
555 : poDS->nStartOfData
556 : + poDS->nRasterXSize*poDS->nRasterYSize*8*iBand,
557 : 8, poDS->nRasterXSize * 8,
558 : GDT_Float64, !CPL_IS_LSB,
559 5 : TRUE, FALSE );
560 :
561 10 : poDS->SetBand( iBand+1, poBand );
562 : }
563 :
564 : /* -------------------------------------------------------------------- */
565 : /* Initialize any PAM information. */
566 : /* -------------------------------------------------------------------- */
567 8 : poDS->SetDescription( poOpenInfo->pszFilename );
568 8 : poDS->TryLoadXML();
569 :
570 : /* -------------------------------------------------------------------- */
571 : /* Check for overviews. */
572 : /* -------------------------------------------------------------------- */
573 8 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
574 :
575 8 : return( poDS );
576 : }
577 :
578 : /************************************************************************/
579 : /* GDALRegister_R() */
580 : /************************************************************************/
581 :
582 338 : void GDALRegister_R()
583 :
584 : {
585 : GDALDriver *poDriver;
586 :
587 338 : if( GDALGetDriverByName( "R" ) == NULL )
588 : {
589 336 : poDriver = new GDALDriver();
590 :
591 336 : poDriver->SetDescription( "R" );
592 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
593 336 : "R Object Data Store" );
594 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
595 336 : "frmt_r.html" );
596 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "rda" );
597 336 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Float32" );
598 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
599 : "<CreationOptionList>"
600 : " <Option name='ASCII' type='boolean' description='For ASCII output, default NO'/>"
601 : " <Option name='COMPRESS' type='boolean' description='Produced Compressed output, default YES'/>"
602 336 : "</CreationOptionList>" );
603 :
604 336 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
605 :
606 336 : poDriver->pfnOpen = RDataset::Open;
607 336 : poDriver->pfnIdentify = RDataset::Identify;
608 336 : poDriver->pfnCreateCopy = RCreateCopy;
609 :
610 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
611 : }
612 338 : }
613 :
|
