1 : /******************************************************************************
2 : *
3 : * Project: ILWIS Driver
4 : * Purpose: GDALDataset driver for ILWIS translator for read/write support.
5 : * Author: Lichun Wang, lichun@itc.nl
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, ITC
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 :
30 : #include "ilwisdataset.h"
31 : #include <float.h>
32 : #include <limits.h>
33 :
34 : // IniFile.cpp: implementation of the IniFile class.
35 : //
36 : //////////////////////////////////////////////////////////////////////
37 0 : bool CompareAsNum::operator() (const string& s1, const string& s2) const
38 : {
39 0 : long Num1 = atoi(s1.c_str());
40 0 : long Num2 = atoi(s2.c_str());
41 0 : return Num1 < Num2;
42 : }
43 :
44 4811 : static string TrimSpaces(const string& input)
45 : {
46 : // find first non space
47 4811 : if ( input.empty())
48 0 : return string();
49 :
50 4811 : size_t iFirstNonSpace = input.find_first_not_of(' ');
51 4811 : size_t iFindLastSpace = input.find_last_not_of(' ');
52 4811 : if (iFirstNonSpace == string::npos || iFindLastSpace == string::npos)
53 0 : return string();
54 :
55 4811 : return input.substr(iFirstNonSpace, iFindLastSpace - iFirstNonSpace + 1);
56 : }
57 :
58 14528 : static string GetLine(FILE* fil)
59 : {
60 14528 : const char *p = CPLReadLineL( fil );
61 14528 : if (p == NULL)
62 34 : return string();
63 :
64 14494 : CPLString osWrk = p;
65 14494 : osWrk.Trim();
66 14494 : return string(osWrk);
67 : }
68 :
69 : //////////////////////////////////////////////////////////////////////
70 : // Construction/Destruction
71 : //////////////////////////////////////////////////////////////////////
72 1195 : IniFile::IniFile(const string& filenam)
73 : {
74 1195 : filename = filenam;
75 1195 : Load();
76 1195 : bChanged = false; // Start tracking changes
77 1195 : }
78 :
79 1195 : IniFile::~IniFile()
80 : {
81 1195 : if (bChanged)
82 : {
83 859 : Store();
84 859 : bChanged = false;
85 : }
86 :
87 4961 : for (Sections::iterator iter = sections.begin(); iter != sections.end(); ++iter)
88 : {
89 3766 : (*(*iter).second).clear();
90 3766 : delete (*iter).second;
91 : }
92 :
93 1195 : sections.clear();
94 1195 : }
95 :
96 8447 : void IniFile::SetKeyValue(const string& section, const string& key, const string& value)
97 : {
98 8447 : Sections::iterator iterSect = sections.find(section);
99 8447 : if (iterSect == sections.end())
100 : {
101 : // Add a new section, with one new key/value entry
102 3766 : SectionEntries *entries = new SectionEntries;
103 3766 : (*entries)[key] = value;
104 3766 : sections[section] = entries;
105 : }
106 : else
107 : {
108 : // Add one new key/value entry in an existing section
109 4681 : SectionEntries *entries = (*iterSect).second;
110 4681 : (*entries)[key] = value;
111 : }
112 8447 : bChanged = true;
113 8447 : }
114 :
115 336 : string IniFile::GetKeyValue(const string& section, const string& key)
116 : {
117 336 : Sections::iterator iterSect = sections.find(section);
118 336 : if (iterSect != sections.end())
119 : {
120 323 : SectionEntries *entries = (*iterSect).second;
121 323 : SectionEntries::iterator iterEntry = (*entries).find(key);
122 323 : if (iterEntry != (*entries).end())
123 276 : return (*iterEntry).second;
124 : }
125 :
126 60 : return string();
127 : }
128 :
129 0 : void IniFile::RemoveKeyValue(const string& section, const string& key)
130 : {
131 0 : Sections::iterator iterSect = sections.find(section);
132 0 : if (iterSect != sections.end())
133 : {
134 : // The section exists, now erase entry "key"
135 0 : SectionEntries *entries = (*iterSect).second;
136 0 : (*entries).erase(key);
137 0 : bChanged = true;
138 : }
139 0 : }
140 :
141 0 : void IniFile::RemoveSection(const string& section)
142 : {
143 0 : Sections::iterator iterSect = sections.find(section);
144 0 : if (iterSect != sections.end())
145 : {
146 : // The section exists, so remove it and all its entries.
147 0 : SectionEntries *entries = (*iterSect).second;
148 0 : (*entries).clear();
149 0 : sections.erase(iterSect);
150 0 : bChanged = true;
151 : }
152 0 : }
153 :
154 1195 : void IniFile::Load()
155 : {
156 : enum ParseState { FindSection, FindKey, ReadFindKey, StoreKey, None } state;
157 1195 : FILE *filIni = VSIFOpenL(filename.c_str(), "r");
158 1195 : if (filIni == NULL)
159 96 : return;
160 :
161 1099 : string section, key, value;
162 1099 : state = FindSection;
163 1099 : string s;
164 28432 : while (!VSIFEofL(filIni) || !s.empty() )
165 : {
166 26234 : switch (state)
167 : {
168 : case FindSection:
169 11058 : s = GetLine(filIni);
170 11058 : if (s.empty())
171 3470 : continue;
172 :
173 7588 : if (s[0] == '[')
174 : {
175 3470 : size_t iLast = s.find_first_of(']');
176 3470 : if (iLast != string::npos)
177 : {
178 3470 : section = s.substr(1, iLast - 1);
179 3470 : state = ReadFindKey;
180 : }
181 : }
182 : else
183 4118 : state = FindKey;
184 7588 : break;
185 : case ReadFindKey:
186 3470 : s = GetLine(filIni); // fall through (no break)
187 : case FindKey:
188 : {
189 7588 : size_t iEqu = s.find_first_of('=');
190 7588 : if (iEqu != string::npos)
191 : {
192 7588 : key = s.substr(0, iEqu);
193 7588 : value = s.substr(iEqu + 1);
194 7588 : state = StoreKey;
195 : }
196 : else
197 0 : state = ReadFindKey;
198 : }
199 7588 : break;
200 : case StoreKey:
201 7588 : SetKeyValue(section, key, value);
202 7588 : state = FindSection;
203 : break;
204 :
205 : case None:
206 : // Do we need to do anything? Perhaps this never occurs.
207 : break;
208 : }
209 : }
210 :
211 1099 : VSIFCloseL(filIni);
212 : }
213 :
214 859 : void IniFile::Store()
215 : {
216 859 : FILE *filIni = VSIFOpenL(filename.c_str(), "w+");
217 859 : if (filIni == NULL)
218 0 : return;
219 :
220 859 : Sections::iterator iterSect;
221 3428 : for (iterSect = sections.begin(); iterSect != sections.end(); ++iterSect)
222 : {
223 2569 : CPLString osLine;
224 :
225 : // write the section name
226 2569 : osLine.Printf( "[%s]\r\n", (*iterSect).first.c_str());
227 2569 : VSIFWriteL( osLine.c_str(), 1, strlen(osLine), filIni );
228 2569 : SectionEntries *entries = (*iterSect).second;
229 2569 : SectionEntries::iterator iterEntry;
230 7380 : for (iterEntry = (*entries).begin(); iterEntry != (*entries).end(); ++iterEntry)
231 : {
232 4811 : string key = (*iterEntry).first;
233 : osLine.Printf( "%s=%s\r\n",
234 4811 : TrimSpaces(key).c_str(), (*iterEntry).second.c_str());
235 4811 : VSIFWriteL( osLine.c_str(), 1, strlen(osLine), filIni );
236 : }
237 :
238 2569 : VSIFWriteL( "\r\n", 1, 2, filIni );
239 : }
240 :
241 859 : VSIFCloseL( filIni );
242 : }
243 :
244 : // End of the implementation of IniFile class. ///////////////////////
245 : //////////////////////////////////////////////////////////////////////
246 :
247 0 : static long longConv(double x) {
248 0 : if ((x == rUNDEF) || (x > LONG_MAX) || (x < LONG_MIN))
249 0 : return iUNDEF;
250 : else
251 0 : return (long)floor(x + 0.5);
252 : }
253 :
254 336 : string ReadElement(string section, string entry, string filename)
255 : {
256 336 : if (section.length() == 0)
257 0 : return string();
258 336 : if (entry.length() == 0)
259 0 : return string();
260 336 : if (filename.length() == 0)
261 0 : return string();
262 :
263 336 : IniFile MyIniFile (filename);
264 :
265 336 : return MyIniFile.GetKeyValue(section, entry);;
266 : }
267 :
268 859 : bool WriteElement(string sSection, string sEntry,
269 : string fn, string sValue)
270 : {
271 859 : if (0 == fn.length())
272 0 : return false;
273 :
274 859 : IniFile MyIniFile (fn);
275 :
276 859 : MyIniFile.SetKeyValue(sSection, sEntry, sValue);
277 859 : return true;
278 : }
279 :
280 53 : bool WriteElement(string sSection, string sEntry,
281 : string fn, int nValue)
282 : {
283 53 : if (0 == fn.length())
284 0 : return false;
285 :
286 : char strdouble[45];
287 53 : sprintf(strdouble, "%d", nValue);
288 53 : string sValue = string(strdouble);
289 106 : return WriteElement(sSection, sEntry, fn, sValue) != 0;
290 : }
291 :
292 78 : bool WriteElement(string sSection, string sEntry,
293 : string fn, double dValue)
294 : {
295 78 : if (0 == fn.length())
296 0 : return false;
297 :
298 : char strdouble[45];
299 78 : sprintf(strdouble, "%.6f", dValue);
300 78 : string sValue = string(strdouble);
301 156 : return WriteElement(sSection, sEntry, fn, sValue) != 0;
302 : }
303 :
304 10 : static CPLErr GetRowCol(string str,int &Row, int &Col)
305 : {
306 10 : string delimStr = " ,;";
307 20 : size_t iPos = str.find_first_of(delimStr);
308 10 : if (iPos != string::npos)
309 : {
310 10 : Row = atoi(str.substr(0, iPos).c_str());
311 : }
312 : else
313 : {
314 : CPLError( CE_Failure, CPLE_AppDefined,
315 0 : "Read of RowCol failed.");
316 0 : return CE_Failure;
317 : }
318 20 : iPos = str.find_last_of(delimStr);
319 10 : if (iPos != string::npos)
320 : {
321 10 : Col = atoi(str.substr(iPos+1, str.length()-iPos).c_str());
322 : }
323 10 : return CE_None;
324 : }
325 :
326 : //! Converts ILWIS data type to GDAL data type.
327 49 : static GDALDataType ILWIS2GDALType(ilwisStoreType stStoreType)
328 : {
329 49 : GDALDataType eDataType = GDT_Unknown;
330 :
331 49 : switch (stStoreType){
332 : case stByte: {
333 36 : eDataType = GDT_Byte;
334 36 : break;
335 : }
336 : case stInt:{
337 4 : eDataType = GDT_Int16;
338 4 : break;
339 : }
340 : case stLong:{
341 4 : eDataType = GDT_Int32;
342 4 : break;
343 : }
344 : case stFloat:{
345 3 : eDataType = GDT_Float32;
346 3 : break;
347 : }
348 : case stReal:{
349 2 : eDataType = GDT_Float64;
350 : break;
351 : }
352 : default: {
353 : break;
354 : }
355 : }
356 :
357 49 : return eDataType;
358 : }
359 :
360 : //Determine store type of ILWIS raster
361 52 : static string GDALType2ILWIS(GDALDataType type)
362 : {
363 52 : string sStoreType;
364 52 : sStoreType = "";
365 52 : switch( type )
366 : {
367 : case GDT_Byte:{
368 33 : sStoreType = "Byte";
369 33 : break;
370 : }
371 : case GDT_Int16:
372 : case GDT_UInt16:{
373 6 : sStoreType = "Int";
374 6 : break;
375 : }
376 : case GDT_Int32:
377 : case GDT_UInt32:{
378 6 : sStoreType = "Long";
379 6 : break;
380 : }
381 : case GDT_Float32:{
382 4 : sStoreType = "Float";
383 4 : break;
384 : }
385 : case GDT_Float64:{
386 3 : sStoreType = "Real";
387 3 : break;
388 : }
389 : default:{
390 : CPLError( CE_Failure, CPLE_NotSupported,
391 : "Data type %s not supported by ILWIS format.\n",
392 0 : GDALGetDataTypeName( type ) );
393 : break;
394 : }
395 : }
396 0 : return sStoreType;
397 : }
398 :
399 68 : static CPLErr GetStoreType(string pszFileName, ilwisStoreType &stStoreType)
400 : {
401 68 : string st = ReadElement("MapStore", "Type", pszFileName.c_str());
402 :
403 136 : if( EQUAL(st.c_str(),"byte"))
404 : {
405 51 : stStoreType = stByte;
406 : }
407 17 : else if( EQUAL(st.c_str(),"int"))
408 : {
409 4 : stStoreType = stInt;
410 : }
411 13 : else if( EQUAL(st.c_str(),"long"))
412 : {
413 4 : stStoreType = stLong;
414 : }
415 9 : else if( EQUAL(st.c_str(),"float"))
416 : {
417 7 : stStoreType = stFloat;
418 : }
419 2 : else if( EQUAL(st.c_str(),"real"))
420 : {
421 2 : stStoreType = stReal;
422 : }
423 : else
424 : {
425 : CPLError( CE_Failure, CPLE_AppDefined,
426 0 : "Unsupported ILWIS store type.");
427 0 : return CE_Failure;
428 : }
429 68 : return CE_None;
430 : }
431 :
432 :
433 39 : ILWISDataset::ILWISDataset()
434 :
435 : {
436 39 : bGeoDirty = FALSE;
437 39 : bNewDataset = FALSE;
438 39 : pszProjection = CPLStrdup("");
439 39 : adfGeoTransform[0] = 0.0;
440 39 : adfGeoTransform[1] = 1.0;
441 39 : adfGeoTransform[2] = 0.0;
442 39 : adfGeoTransform[3] = 0.0;
443 39 : adfGeoTransform[4] = 0.0;
444 39 : adfGeoTransform[5] = 1.0;
445 39 : }
446 :
447 : /************************************************************************/
448 : /* ~ILWISDataset() */
449 : /************************************************************************/
450 :
451 78 : ILWISDataset::~ILWISDataset()
452 :
453 : {
454 39 : FlushCache();
455 39 : CPLFree( pszProjection );
456 78 : }
457 :
458 : /************************************************************************/
459 : /* CollectTransformCoef() */
460 : /* */
461 : /* Collect the geotransform, support for the GeoRefCorners */
462 : /* georeferencing only; We use the extent of the coordinates */
463 : /* to determine the pixelsize in X and Y direction. Then calculate */
464 : /* the transform coefficients from the extent and pixelsize */
465 : /************************************************************************/
466 :
467 10 : void ILWISDataset::CollectTransformCoef(string &pszRefName)
468 :
469 : {
470 10 : pszRefName = "";
471 10 : string georef;
472 10 : if ( EQUAL(pszFileType.c_str(),"Map") )
473 7 : georef = ReadElement("Map", "GeoRef", osFileName);
474 : else
475 3 : georef = ReadElement("MapList", "GeoRef", osFileName);
476 :
477 : //Capture the geotransform, only if the georef is not 'none',
478 : //otherwise, the default transform should be returned.
479 10 : if( (georef.length() != 0) && !EQUAL(georef.c_str(),"none"))
480 : {
481 : //Form the geo-referencing name
482 10 : string pszBaseName = string(CPLGetBasename(georef.c_str()) );
483 20 : string pszPath = string(CPLGetPath( osFileName ));
484 : pszRefName = string(CPLFormFilename(pszPath.c_str(),
485 20 : pszBaseName.c_str(),"grf" ));
486 :
487 : //Check the geo-reference type,support for the GeoRefCorners only
488 20 : string georeftype = ReadElement("GeoRef", "Type", pszRefName);
489 10 : if (EQUAL(georeftype.c_str(),"GeoRefCorners"))
490 : {
491 : //Center or top-left corner of the pixel approach?
492 5 : string IsCorner = ReadElement("GeoRefCorners", "CornersOfCorners", pszRefName);
493 :
494 : //Collect the extent of the coordinates
495 5 : string sMinX = ReadElement("GeoRefCorners", "MinX", pszRefName);
496 5 : string sMinY = ReadElement("GeoRefCorners", "MinY", pszRefName);
497 5 : string sMaxX = ReadElement("GeoRefCorners", "MaxX", pszRefName);
498 5 : string sMaxY = ReadElement("GeoRefCorners", "MaxY", pszRefName);
499 :
500 : //Calculate pixel size in X and Y direction from the extent
501 5 : double deltaX = atof(sMaxX.c_str()) - atof(sMinX.c_str());
502 5 : double deltaY = atof(sMaxY.c_str()) - atof(sMinY.c_str());
503 :
504 5 : double PixelSizeX = deltaX / (double)nRasterXSize;
505 5 : double PixelSizeY = deltaY / (double)nRasterYSize;
506 :
507 5 : if (EQUAL(IsCorner.c_str(),"Yes"))
508 : {
509 5 : adfGeoTransform[0] = atof(sMinX.c_str());
510 5 : adfGeoTransform[3] = atof(sMaxY.c_str());
511 : }
512 : else
513 : {
514 0 : adfGeoTransform[0] = atof(sMinX.c_str()) - PixelSizeX/2.0;
515 0 : adfGeoTransform[3] = atof(sMaxY.c_str()) + PixelSizeY/2.0;
516 : }
517 :
518 5 : adfGeoTransform[1] = PixelSizeX;
519 5 : adfGeoTransform[2] = 0.0;
520 5 : adfGeoTransform[4] = 0.0;
521 5 : adfGeoTransform[5] = -PixelSizeY;
522 10 : }
523 :
524 10 : }
525 :
526 10 : }
527 :
528 : /************************************************************************/
529 : /* WriteGeoReference() */
530 : /* */
531 : /* Try to write a geo-reference file for the dataset to create */
532 : /************************************************************************/
533 :
534 13 : CPLErr ILWISDataset::WriteGeoReference()
535 : {
536 : //check wheather we should write out a georeference file.
537 : //dataset must be north up
538 13 : if( adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
539 0 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
540 0 : || adfGeoTransform[4] != 0.0 || fabs(adfGeoTransform[5]) != 1.0 )
541 : {
542 13 : SetGeoTransform( adfGeoTransform ); // is this needed?
543 13 : if (adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0)
544 : {
545 13 : int nXSize = GetRasterXSize();
546 13 : int nYSize = GetRasterYSize();
547 13 : double dLLLat = (adfGeoTransform[3]
548 13 : + nYSize * adfGeoTransform[5] );
549 13 : double dLLLong = (adfGeoTransform[0] );
550 13 : double dURLat = (adfGeoTransform[3] );
551 13 : double dURLong = (adfGeoTransform[0]
552 13 : + nXSize * adfGeoTransform[1] );
553 :
554 13 : string grFileName = CPLResetExtension(osFileName, "grf" );
555 26 : WriteElement("Ilwis", "Type", grFileName, "GeoRef");
556 26 : WriteElement("GeoRef", "lines", grFileName, nYSize);
557 13 : WriteElement("GeoRef", "columns", grFileName, nXSize);
558 13 : WriteElement("GeoRef", "Type", grFileName, "GeoRefCorners");
559 26 : WriteElement("GeoRefCorners", "CornersOfCorners", grFileName, "Yes");
560 26 : WriteElement("GeoRefCorners", "MinX", grFileName, dLLLong);
561 13 : WriteElement("GeoRefCorners", "MinY", grFileName, dLLLat);
562 13 : WriteElement("GeoRefCorners", "MaxX", grFileName, dURLong);
563 13 : WriteElement("GeoRefCorners", "MaxY", grFileName, dURLat);
564 :
565 : //Re-write the GeoRef property to raster ODF
566 : //Form band file name
567 13 : string sBaseName = string(CPLGetBasename(osFileName) );
568 26 : string sPath = string(CPLGetPath(osFileName));
569 13 : if (nBands == 1)
570 : {
571 9 : WriteElement("Map", "GeoRef", osFileName, sBaseName + ".grf");
572 : }
573 : else
574 : {
575 18 : for( int iBand = 0; iBand < nBands; iBand++ )
576 : {
577 14 : if (iBand == 0)
578 4 : WriteElement("MapList", "GeoRef", osFileName, sBaseName + ".grf");
579 : char pszName[100];
580 14 : sprintf(pszName, "%s_band_%d", sBaseName.c_str(),iBand + 1 );
581 14 : string pszODFName = string(CPLFormFilename(sPath.c_str(),pszName,"mpr"));
582 28 : WriteElement("Map", "GeoRef", pszODFName, sBaseName + ".grf");
583 : }
584 13 : }
585 : }
586 : }
587 13 : return CE_None;
588 : }
589 :
590 : /************************************************************************/
591 : /* GetProjectionRef() */
592 : /************************************************************************/
593 :
594 32 : const char *ILWISDataset::GetProjectionRef()
595 :
596 : {
597 32 : return ( pszProjection );
598 : }
599 :
600 : /************************************************************************/
601 : /* SetProjection() */
602 : /************************************************************************/
603 :
604 13 : CPLErr ILWISDataset::SetProjection( const char * pszNewProjection )
605 :
606 : {
607 13 : CPLFree( pszProjection );
608 13 : pszProjection = CPLStrdup( pszNewProjection );
609 13 : bGeoDirty = TRUE;
610 :
611 13 : return CE_None;
612 : }
613 :
614 : /************************************************************************/
615 : /* GetGeoTransform() */
616 : /************************************************************************/
617 :
618 20 : CPLErr ILWISDataset::GetGeoTransform( double * padfTransform )
619 :
620 : {
621 20 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
622 20 : return( CE_None );
623 : }
624 :
625 : /************************************************************************/
626 : /* SetGeoTransform() */
627 : /************************************************************************/
628 :
629 26 : CPLErr ILWISDataset::SetGeoTransform( double * padfTransform )
630 :
631 : {
632 26 : memmove( adfGeoTransform, padfTransform, sizeof(double)*6 );
633 :
634 26 : if (adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0)
635 26 : bGeoDirty = TRUE;
636 :
637 26 : return CE_None;
638 : }
639 :
640 10 : bool CheckASCII(unsigned char * buf, int size)
641 : {
642 2539 : for (int i = 0; i < size; ++i)
643 2529 : if (!isascii(buf[i]))
644 0 : return false;
645 :
646 10 : return true;
647 : }
648 : /************************************************************************/
649 : /* Open() */
650 : /************************************************************************/
651 :
652 9538 : GDALDataset *ILWISDataset::Open( GDALOpenInfo * poOpenInfo )
653 :
654 : {
655 : /* -------------------------------------------------------------------- */
656 : /* Does this look like an ILWIS file */
657 : /* -------------------------------------------------------------------- */
658 9538 : if( poOpenInfo->nHeaderBytes < 1 )
659 8529 : return NULL;
660 :
661 1009 : string sExt = CPLGetExtension( poOpenInfo->pszFilename );
662 2018 : if (!EQUAL(sExt.c_str(),"mpr") && !EQUAL(sExt.c_str(),"mpl"))
663 999 : return NULL;
664 :
665 10 : if (!CheckASCII(poOpenInfo->pabyHeader, poOpenInfo->nHeaderBytes))
666 0 : return NULL;
667 :
668 10 : string ilwistype = ReadElement("Ilwis", "Type", poOpenInfo->pszFilename);
669 20 : if( ilwistype.length() == 0)
670 0 : return NULL;
671 :
672 10 : string sFileType; //map or map list
673 : int iBandCount;
674 10 : string mapsize;
675 10 : string maptype = ReadElement("BaseMap", "Type", poOpenInfo->pszFilename);
676 20 : string sBaseName = string(CPLGetBasename(poOpenInfo->pszFilename) );
677 20 : string sPath = string(CPLGetPath( poOpenInfo->pszFilename));
678 :
679 : //Verify whether it is a map list or a map
680 20 : if( EQUAL(ilwistype.c_str(),"MapList") )
681 : {
682 3 : sFileType = string("MapList");
683 6 : string sMaps = ReadElement("MapList", "Maps", poOpenInfo->pszFilename);
684 6 : iBandCount = atoi(sMaps.c_str());
685 3 : mapsize = ReadElement("MapList", "Size", poOpenInfo->pszFilename);
686 3 : for (int iBand = 0; iBand < iBandCount; ++iBand )
687 : {
688 : //Form the band file name.
689 : char cBandName[45];
690 6 : sprintf( cBandName, "Map%d", iBand);
691 6 : string sBandName = ReadElement("MapList", string(cBandName), poOpenInfo->pszFilename);
692 12 : string pszBandBaseName = string(CPLGetBasename(sBandName.c_str()) );
693 12 : string pszBandPath = string(CPLGetPath( sBandName.c_str()));
694 12 : if ( 0 == pszBandPath.length() )
695 : {
696 : sBandName = string(CPLFormFilename(sPath.c_str(),
697 6 : pszBandBaseName.c_str(),"mpr" ));
698 : }
699 : //Verify the file exetension, it must be an ILWIS raw data file
700 : //with extension .mp#, otherwise, unsupported
701 : //This drive only supports a map list which stores a set of ILWIS raster maps,
702 6 : string sMapStoreName = ReadElement("MapStore", "Data", sBandName);
703 6 : string sExt = CPLGetExtension( sMapStoreName.c_str() );
704 12 : if ( !EQUALN( sExt.c_str(), "mp#", 3 ))
705 : {
706 : CPLError( CE_Failure, CPLE_AppDefined,
707 : "Unsupported ILWIS data file. \n"
708 0 : "can't treat as raster.\n" );
709 0 : return FALSE;
710 : }
711 0 : }
712 : }
713 7 : else if(EQUAL(ilwistype.c_str(),"BaseMap") && EQUAL(maptype.c_str(),"Map"))
714 : {
715 7 : sFileType = "Map";
716 7 : iBandCount = 1;
717 7 : mapsize = ReadElement("Map", "Size", poOpenInfo->pszFilename);
718 14 : string sMapType = ReadElement("Map", "Type", poOpenInfo->pszFilename);
719 : ilwisStoreType stStoreType;
720 14 : if (
721 : GetStoreType(string(poOpenInfo->pszFilename), stStoreType) != CE_None )
722 : {
723 : //CPLError( CE_Failure, CPLE_AppDefined,
724 : // "Unsupported ILWIS data file. \n"
725 : // "can't treat as raster.\n" );
726 0 : return FALSE;
727 0 : }
728 : }
729 : else
730 : {
731 : CPLError( CE_Failure, CPLE_AppDefined,
732 : "Unsupported ILWIS data file. \n"
733 0 : "can't treat as raster.\n" );
734 0 : return FALSE;
735 : }
736 :
737 : /* -------------------------------------------------------------------- */
738 : /* Create a corresponding GDALDataset. */
739 : /* -------------------------------------------------------------------- */
740 : ILWISDataset *poDS;
741 10 : poDS = new ILWISDataset();
742 :
743 : /* -------------------------------------------------------------------- */
744 : /* Capture raster size from ILWIS file (.mpr). */
745 : /* -------------------------------------------------------------------- */
746 10 : int Row = 0, Col = 0;
747 20 : if ( GetRowCol(mapsize, Row, Col) != CE_None)
748 0 : return FALSE;
749 10 : poDS->nRasterXSize = Col;
750 10 : poDS->nRasterYSize = Row;
751 10 : poDS->osFileName = poOpenInfo->pszFilename;
752 10 : poDS->pszFileType = sFileType;
753 : /* -------------------------------------------------------------------- */
754 : /* Create band information objects. */
755 : /* -------------------------------------------------------------------- */
756 : //poDS->pszFileName = new char[strlen(poOpenInfo->pszFilename) + 1];
757 10 : poDS->nBands = iBandCount;
758 23 : for( int iBand = 0; iBand < poDS->nBands; iBand++ )
759 : {
760 13 : poDS->SetBand( iBand+1, new ILWISRasterBand( poDS, iBand+1 ) );
761 : }
762 :
763 : /* -------------------------------------------------------------------- */
764 : /* Collect the geotransform coefficients */
765 : /* -------------------------------------------------------------------- */
766 10 : string pszGeoRef;
767 10 : poDS->CollectTransformCoef(pszGeoRef);
768 :
769 : /* -------------------------------------------------------------------- */
770 : /* Translation from ILWIS coordinate system definition */
771 : /* -------------------------------------------------------------------- */
772 10 : if( (pszGeoRef.length() != 0) && !EQUAL(pszGeoRef.c_str(),"none"))
773 : {
774 :
775 : // Fetch coordinate system
776 10 : string csy = ReadElement("GeoRef", "CoordSystem", pszGeoRef);
777 10 : string pszProj;
778 :
779 10 : if( (csy.length() != 0) && !EQUAL(csy.c_str(),"unknown.csy"))
780 : {
781 :
782 : //Form the coordinate system file name
783 5 : if( !(EQUALN( csy.c_str(), "latlon.csy", 10 )) &&
784 : !(EQUALN( csy.c_str(), "LatlonWGS84.csy", 15 )))
785 : {
786 5 : string pszBaseName = string(CPLGetBasename(csy.c_str()) );
787 10 : string pszPath = string(CPLGetPath( poDS->osFileName ));
788 : csy = string(CPLFormFilename(pszPath.c_str(),
789 10 : pszBaseName.c_str(),"csy" ));
790 10 : pszProj = ReadElement("CoordSystem", "Type", csy);
791 5 : if (pszProj.length() == 0 ) //default to projection
792 0 : pszProj = "Projection";
793 : }
794 : else
795 : {
796 0 : pszProj = "LatLon";
797 : }
798 :
799 5 : if( (EQUALN( pszProj.c_str(), "LatLon", 6 )) ||
800 : (EQUALN( pszProj.c_str(), "Projection", 10 )))
801 5 : poDS->ReadProjection( csy );
802 10 : }
803 : }
804 :
805 : /* -------------------------------------------------------------------- */
806 : /* Initialize any PAM information. */
807 : /* -------------------------------------------------------------------- */
808 10 : poDS->SetDescription( poOpenInfo->pszFilename );
809 10 : poDS->TryLoadXML();
810 :
811 10 : return( poDS );
812 : }
813 :
814 : /************************************************************************/
815 : /* FlushCache() */
816 : /************************************************************************/
817 :
818 53 : void ILWISDataset::FlushCache()
819 :
820 : {
821 53 : GDALDataset::FlushCache();
822 :
823 53 : if( bGeoDirty == TRUE )
824 : {
825 13 : WriteGeoReference();
826 13 : WriteProjection();
827 13 : bGeoDirty = FALSE;
828 : }
829 53 : }
830 :
831 : /************************************************************************/
832 : /* Create() */
833 : /* */
834 : /* Create a new ILWIS file. */
835 : /************************************************************************/
836 :
837 37 : GDALDataset *ILWISDataset::Create(const char* pszFilename,
838 : int nXSize, int nYSize,
839 : int nBands, GDALDataType eType,
840 : char** papszParmList)
841 : {
842 : ILWISDataset *poDS;
843 : int iBand;
844 :
845 : /* -------------------------------------------------------------------- */
846 : /* Verify input options. */
847 : /* -------------------------------------------------------------------- */
848 37 : if( eType != GDT_Byte && eType != GDT_Int16 && eType != GDT_Int32
849 : && eType != GDT_Float32 && eType != GDT_Float64 && eType != GDT_UInt16 && eType != GDT_UInt32)
850 : {
851 : CPLError( CE_Failure, CPLE_AppDefined,
852 : "Attempt to create ILWIS dataset with an illegal\n"
853 : "data type (%s).\n",
854 8 : GDALGetDataTypeName(eType) );
855 :
856 8 : return NULL;
857 : }
858 :
859 : /* -------------------------------------------------------------------- */
860 : /* Translate the data type. */
861 : /* Determine store type of ILWIS raster */
862 : /* -------------------------------------------------------------------- */
863 29 : string sDomain= "value.dom";
864 29 : double stepsize = 1;
865 29 : string sStoreType = GDALType2ILWIS(eType);
866 29 : if( EQUAL(sStoreType.c_str(),""))
867 0 : return NULL;
868 29 : else if( EQUAL(sStoreType.c_str(),"Real") || EQUAL(sStoreType.c_str(),"float"))
869 5 : stepsize = 0;
870 :
871 29 : string pszBaseName = string(CPLGetBasename( pszFilename ));
872 58 : string pszPath = string(CPLGetPath( pszFilename ));
873 :
874 : /* -------------------------------------------------------------------- */
875 : /* Write out object definition file for each band */
876 : /* -------------------------------------------------------------------- */
877 29 : string pszODFName;
878 29 : string pszDataBaseName;
879 29 : string pszFileName;
880 :
881 : char strsize[45];
882 29 : sprintf(strsize, "%d %d", nYSize, nXSize);
883 :
884 : //Form map/maplist name.
885 29 : if ( nBands == 1 )
886 : {
887 17 : pszODFName = string(CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpr"));
888 17 : pszDataBaseName = pszBaseName;
889 17 : pszFileName = CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpr");
890 : }
891 : else
892 : {
893 12 : pszFileName = CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpl");
894 12 : WriteElement("Ilwis", "Type", string(pszFileName), "MapList");
895 24 : WriteElement("MapList", "GeoRef", string(pszFileName), "none.grf");
896 24 : WriteElement("MapList", "Size", string(pszFileName), string(strsize));
897 24 : WriteElement("MapList", "Maps", string(pszFileName), nBands);
898 : }
899 :
900 29 : for( iBand = 0; iBand < nBands; iBand++ )
901 : {
902 48 : if ( nBands > 1 )
903 : {
904 : char pszBandName[100];
905 31 : sprintf(pszBandName, "%s_band_%d", pszBaseName.c_str(),iBand + 1 );
906 31 : pszODFName = string(pszBandName) + ".mpr";
907 62 : pszDataBaseName = string(pszBandName);
908 31 : sprintf(pszBandName, "Map%d", iBand);
909 62 : WriteElement("MapList", string(pszBandName), string(pszFileName), pszODFName);
910 31 : pszODFName = CPLFormFilename(pszPath.c_str(),pszDataBaseName.c_str(),"mpr");
911 : }
912 : /* -------------------------------------------------------------------- */
913 : /* Write data definition per band (.mpr) */
914 : /* -------------------------------------------------------------------- */
915 :
916 48 : WriteElement("Ilwis", "Type", pszODFName, "BaseMap");
917 96 : WriteElement("BaseMap", "Type", pszODFName, "Map");
918 96 : WriteElement("Map", "Type", pszODFName, "MapStore");
919 :
920 96 : WriteElement("BaseMap", "Domain", pszODFName, sDomain);
921 48 : string pszDataName = pszDataBaseName + ".mp#";
922 48 : WriteElement("MapStore", "Data", pszODFName, pszDataName);
923 48 : WriteElement("MapStore", "Structure", pszODFName, "Line");
924 : // sStoreType is used by ILWISRasterBand constructor to determine eDataType
925 96 : WriteElement("MapStore", "Type", pszODFName, sStoreType);
926 :
927 : // For now write-out a "Range" that is as broad as possible.
928 : // If later a better range is found (by inspecting metadata in the source dataset),
929 : // the "Range" will be overwritten by a better version.
930 : double adfMinMax[2];
931 48 : adfMinMax[0] = -9999999.9;
932 48 : adfMinMax[1] = 9999999.9;
933 : char strdouble[45];
934 48 : sprintf(strdouble, "%.3f:%.3f:%3f:offset=0", adfMinMax[0], adfMinMax[1],stepsize);
935 48 : string range = string(strdouble);
936 96 : WriteElement("BaseMap", "Range", pszODFName, range);
937 :
938 48 : WriteElement("Map", "GeoRef", pszODFName, "none.grf");
939 96 : WriteElement("Map", "Size", pszODFName, string(strsize));
940 :
941 : /* -------------------------------------------------------------------- */
942 : /* Try to create the data file. */
943 : /* -------------------------------------------------------------------- */
944 96 : pszDataName = CPLResetExtension(pszODFName.c_str(), "mp#" );
945 :
946 48 : FILE *fp = VSIFOpenL( pszDataName.c_str(), "wb" );
947 :
948 48 : if( fp == NULL )
949 : {
950 : CPLError( CE_Failure, CPLE_OpenFailed,
951 0 : "Unable to create file %s.\n", pszDataName.c_str() );
952 0 : return NULL;
953 : }
954 48 : VSIFCloseL( fp );
955 : }
956 29 : poDS = new ILWISDataset();
957 29 : poDS->nRasterXSize = nXSize;
958 29 : poDS->nRasterYSize = nYSize;
959 29 : poDS->nBands = nBands;
960 29 : poDS->eAccess = GA_Update;
961 29 : poDS->bNewDataset = TRUE;
962 29 : poDS->SetDescription(pszFilename);
963 29 : poDS->osFileName = pszFileName;
964 29 : poDS->pszIlwFileName = string(pszFileName);
965 29 : if ( nBands == 1 )
966 17 : poDS->pszFileType = "Map";
967 : else
968 12 : poDS->pszFileType = "MapList";
969 :
970 : /* -------------------------------------------------------------------- */
971 : /* Create band information objects. */
972 : /* -------------------------------------------------------------------- */
973 :
974 77 : for( iBand = 1; iBand <= poDS->nBands; iBand++ )
975 : {
976 48 : poDS->SetBand( iBand, new ILWISRasterBand( poDS, iBand ) );
977 : }
978 :
979 29 : return poDS;
980 : //return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
981 : }
982 :
983 : /************************************************************************/
984 : /* CreateCopy() */
985 : /************************************************************************/
986 :
987 : GDALDataset *
988 18 : ILWISDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
989 : int bStrict, char ** papszOptions,
990 : GDALProgressFunc pfnProgress, void * pProgressData )
991 :
992 : {
993 :
994 : ILWISDataset *poDS;
995 18 : GDALDataType eType = GDT_Byte;
996 : int iBand;
997 : (void) bStrict;
998 :
999 :
1000 18 : int nXSize = poSrcDS->GetRasterXSize();
1001 18 : int nYSize = poSrcDS->GetRasterYSize();
1002 18 : int nBands = poSrcDS->GetRasterCount();
1003 :
1004 18 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1005 0 : return NULL;
1006 :
1007 : /* -------------------------------------------------------------------- */
1008 : /* Create the basic dataset. */
1009 : /* -------------------------------------------------------------------- */
1010 45 : for( iBand = 0; iBand < nBands; iBand++ )
1011 : {
1012 27 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1013 27 : eType = GDALDataTypeUnion( eType, poBand->GetRasterDataType() );
1014 : }
1015 :
1016 : poDS = (ILWISDataset *) Create( pszFilename,
1017 : poSrcDS->GetRasterXSize(),
1018 : poSrcDS->GetRasterYSize(),
1019 : nBands,
1020 18 : eType, papszOptions );
1021 :
1022 18 : if( poDS == NULL )
1023 4 : return NULL;
1024 14 : string pszBaseName = string(CPLGetBasename( pszFilename ));
1025 28 : string pszPath = string(CPLGetPath( pszFilename ));
1026 :
1027 : /* -------------------------------------------------------------------- */
1028 : /* Copy and geo-transform and projection information. */
1029 : /* -------------------------------------------------------------------- */
1030 : double adfGeoTransform[6];
1031 28 : string georef = "none.grf";
1032 :
1033 : //check wheather we should create georeference file.
1034 : //source dataset must be north up
1035 47 : if( poSrcDS->GetGeoTransform( adfGeoTransform ) == CE_None
1036 15 : && (adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
1037 2 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
1038 2 : || adfGeoTransform[4] != 0.0 || fabs(adfGeoTransform[5]) != 1.0))
1039 : {
1040 13 : poDS->SetGeoTransform( adfGeoTransform );
1041 13 : if (adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0)
1042 13 : georef = pszBaseName + ".grf";
1043 : }
1044 :
1045 14 : const char *pszProj = poSrcDS->GetProjectionRef();
1046 14 : if( pszProj != NULL && strlen(pszProj) > 0 )
1047 13 : poDS->SetProjection( pszProj );
1048 :
1049 : /* -------------------------------------------------------------------- */
1050 : /* Create the output raster files for each band */
1051 : /* -------------------------------------------------------------------- */
1052 :
1053 14 : for( iBand = 0; iBand < nBands; iBand++ )
1054 : {
1055 : FILE *fpData;
1056 : GByte *pData;
1057 :
1058 23 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1059 23 : ILWISRasterBand *desBand = (ILWISRasterBand *) poDS->GetRasterBand( iBand+1 );
1060 :
1061 : /* -------------------------------------------------------------------- */
1062 : /* Translate the data type. */
1063 : /* -------------------------------------------------------------------- */
1064 23 : int nLineSize = nXSize * GDALGetDataTypeSize(eType) / 8;
1065 23 : pData = (GByte *) CPLMalloc( nLineSize );
1066 :
1067 : //Determine the nodata value
1068 : int bHasNoDataValue;
1069 23 : double dNoDataValue = poBand->GetNoDataValue(&bHasNoDataValue);
1070 :
1071 : //Determine store type of ILWIS raster
1072 23 : string sStoreType = GDALType2ILWIS( eType );
1073 23 : double stepsize = 1;
1074 23 : if( EQUAL(sStoreType.c_str(),""))
1075 0 : return NULL;
1076 23 : else if( EQUAL(sStoreType.c_str(),"Real") || EQUAL(sStoreType.c_str(),"float"))
1077 2 : stepsize = 0;
1078 :
1079 : //Form the image file name, create the object definition file.
1080 23 : string pszODFName;
1081 23 : string pszDataBaseName;
1082 23 : if (nBands == 1)
1083 : {
1084 9 : pszODFName = string(CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpr"));
1085 9 : pszDataBaseName = pszBaseName;
1086 : }
1087 : else
1088 : {
1089 : char pszName[100];
1090 14 : sprintf(pszName, "%s_band_%d", pszBaseName.c_str(),iBand + 1 );
1091 14 : pszODFName = string(CPLFormFilename(pszPath.c_str(),pszName,"mpr"));
1092 28 : pszDataBaseName = string(pszName);
1093 : }
1094 : /* -------------------------------------------------------------------- */
1095 : /* Write data definition file for each band (.mpr) */
1096 : /* -------------------------------------------------------------------- */
1097 :
1098 : double adfMinMax[2];
1099 : int bGotMin, bGotMax;
1100 :
1101 23 : adfMinMax[0] = poBand->GetMinimum( &bGotMin );
1102 23 : adfMinMax[1] = poBand->GetMaximum( &bGotMax );
1103 23 : if( ! (bGotMin && bGotMax) )
1104 23 : GDALComputeRasterMinMax((GDALRasterBandH)poBand, FALSE, adfMinMax);
1105 23 : if ((!CPLIsNan(adfMinMax[0])) && CPLIsFinite(adfMinMax[0]) && (!CPLIsNan(adfMinMax[1])) && CPLIsFinite(adfMinMax[1]))
1106 : {
1107 : // only write a range if we got a correct one from the source dataset (otherwise ILWIS can't show the map properly)
1108 : char strdouble[45];
1109 23 : sprintf(strdouble, "%.3f:%.3f:%3f:offset=0", adfMinMax[0], adfMinMax[1],stepsize);
1110 23 : string range = string(strdouble);
1111 46 : WriteElement("BaseMap", "Range", pszODFName, range);
1112 : }
1113 23 : WriteElement("Map", "GeoRef", pszODFName, georef);
1114 :
1115 : /* -------------------------------------------------------------------- */
1116 : /* Loop over image, copy the image data. */
1117 : /* -------------------------------------------------------------------- */
1118 23 : CPLErr eErr = CE_None;
1119 :
1120 : //For file name for raw data, and create binary files.
1121 23 : string pszDataFileName = CPLResetExtension(pszODFName.c_str(), "mp#" );
1122 :
1123 23 : fpData = desBand->fpRaw;
1124 23 : if( fpData == NULL )
1125 : {
1126 : CPLError( CE_Failure, CPLE_OpenFailed,
1127 : "Attempt to create file `%s' failed.\n",
1128 0 : pszFilename );
1129 0 : return NULL;
1130 : }
1131 :
1132 273 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ )
1133 : {
1134 : eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
1135 : pData, nXSize, 1, eType,
1136 250 : 0, 0 );
1137 :
1138 250 : if( eErr == CE_None )
1139 : {
1140 250 : if (bHasNoDataValue)
1141 : {
1142 : // pData may have entries with value = dNoDataValue
1143 : // ILWIS uses a fixed value for nodata, depending on the data-type
1144 : // Therefore translate the NoDataValue from each band to ILWIS
1145 0 : for (int iCol = 0; iCol < nXSize; iCol++ )
1146 : {
1147 0 : if( EQUAL(sStoreType.c_str(),"Byte"))
1148 : {
1149 0 : if ( ((GByte * )pData)[iCol] == dNoDataValue )
1150 0 : (( GByte * )pData)[iCol] = 0;
1151 : }
1152 0 : else if( EQUAL(sStoreType.c_str(),"Int"))
1153 : {
1154 0 : if ( ((GInt16 * )pData)[iCol] == dNoDataValue )
1155 0 : (( GInt16 * )pData)[iCol] = shUNDEF;
1156 : }
1157 0 : else if( EQUAL(sStoreType.c_str(),"Long"))
1158 : {
1159 0 : if ( ((GInt32 * )pData)[iCol] == dNoDataValue )
1160 0 : (( GInt32 * )pData)[iCol] = iUNDEF;
1161 : }
1162 0 : else if( EQUAL(sStoreType.c_str(),"float"))
1163 : {
1164 0 : if ((((float * )pData)[iCol] == dNoDataValue ) || (CPLIsNan((( float * )pData)[iCol])))
1165 0 : (( float * )pData)[iCol] = flUNDEF;
1166 : }
1167 0 : else if( EQUAL(sStoreType.c_str(),"Real"))
1168 : {
1169 0 : if ((((double * )pData)[iCol] == dNoDataValue ) || (CPLIsNan((( double * )pData)[iCol])))
1170 0 : (( double * )pData)[iCol] = rUNDEF;
1171 : }
1172 : }
1173 : }
1174 250 : int iSize = VSIFWriteL( pData, 1, nLineSize, desBand->fpRaw );
1175 250 : if ( iSize < 1 )
1176 : {
1177 0 : CPLFree( pData );
1178 : //CPLFree( pData32 );
1179 : CPLError( CE_Failure, CPLE_FileIO,
1180 0 : "Write of file failed with fwrite error.");
1181 0 : return NULL;
1182 : }
1183 : }
1184 250 : if( !pfnProgress(iLine / (nYSize * nBands), NULL, pProgressData ) )
1185 0 : return NULL;
1186 : }
1187 23 : VSIFFlushL( fpData );
1188 23 : CPLFree( pData );
1189 : }
1190 :
1191 14 : poDS->FlushCache();
1192 :
1193 14 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1194 : {
1195 : CPLError( CE_Failure, CPLE_UserInterrupt,
1196 0 : "User terminated" );
1197 0 : delete poDS;
1198 :
1199 : GDALDriver *poILWISDriver =
1200 0 : (GDALDriver *) GDALGetDriverByName( "ILWIS" );
1201 0 : poILWISDriver->Delete( pszFilename );
1202 0 : return NULL;
1203 : }
1204 :
1205 14 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1206 :
1207 14 : return poDS;
1208 : }
1209 :
1210 : /************************************************************************/
1211 : /* ILWISRasterBand() */
1212 : /************************************************************************/
1213 :
1214 61 : ILWISRasterBand::ILWISRasterBand( ILWISDataset *poDS, int nBand )
1215 :
1216 : {
1217 61 : string sBandName;
1218 61 : if ( EQUAL(poDS->pszFileType.c_str(),"Map"))
1219 24 : sBandName = string(poDS->osFileName);
1220 : else //map list
1221 : {
1222 : //Form the band name
1223 : char cBandName[45];
1224 37 : sprintf( cBandName, "Map%d", nBand-1);
1225 37 : sBandName = ReadElement("MapList", string(cBandName), string(poDS->osFileName));
1226 37 : string sInputPath = string(CPLGetPath( poDS->osFileName));
1227 74 : string sBandPath = string(CPLGetPath( sBandName.c_str()));
1228 74 : string sBandBaseName = string(CPLGetBasename( sBandName.c_str()));
1229 74 : if ( 0==sBandPath.length() )
1230 37 : sBandName = string(CPLFormFilename(sInputPath.c_str(),sBandBaseName.c_str(),"mpr" ));
1231 : else
1232 0 : sBandName = string(CPLFormFilename(sBandPath.c_str(),sBandBaseName.c_str(),"mpr" ));
1233 : }
1234 :
1235 61 : if (poDS->bNewDataset)
1236 : {
1237 : // Called from Create():
1238 : // eDataType is defaulted to GDT_Byte by GDALRasterBand::GDALRasterBand
1239 : // Here we set it to match the value of sStoreType (that was set in ILWISDataset::Create)
1240 : // Unfortunately we can't take advantage of the ILWIS "ValueRange" object that would use
1241 : // the most compact storeType possible, without going through all values.
1242 48 : GetStoreType(sBandName, psInfo.stStoreType);
1243 48 : eDataType = ILWIS2GDALType(psInfo.stStoreType);
1244 : }
1245 : else // Called from Open(), thus convert ILWIS type from ODF to eDataType
1246 13 : GetILWISInfo(sBandName);
1247 61 : this->poDS = poDS;
1248 61 : this->nBand = nBand;
1249 61 : nBlockXSize = poDS->GetRasterXSize();
1250 61 : nBlockYSize = 1;
1251 61 : switch (psInfo.stStoreType)
1252 : {
1253 : case stByte:
1254 46 : nSizePerPixel = GDALGetDataTypeSize(GDT_Byte) / 8;
1255 46 : break;
1256 : case stInt:
1257 4 : nSizePerPixel = GDALGetDataTypeSize(GDT_Int16) / 8;
1258 4 : break;
1259 : case stLong:
1260 4 : nSizePerPixel = GDALGetDataTypeSize(GDT_Int32) / 8;
1261 4 : break;
1262 : case stFloat:
1263 5 : nSizePerPixel = GDALGetDataTypeSize(GDT_Float32) / 8;
1264 5 : break;
1265 : case stReal:
1266 2 : nSizePerPixel = GDALGetDataTypeSize(GDT_Float64) / 8;
1267 : break;
1268 : }
1269 61 : ILWISOpen(sBandName);
1270 61 : }
1271 :
1272 : /************************************************************************/
1273 : /* ~ILWISRasterBand() */
1274 : /************************************************************************/
1275 :
1276 122 : ILWISRasterBand::~ILWISRasterBand()
1277 :
1278 : {
1279 61 : if( fpRaw != NULL )
1280 : {
1281 61 : VSIFCloseL( fpRaw );
1282 61 : fpRaw = NULL;
1283 : }
1284 122 : }
1285 :
1286 :
1287 : /************************************************************************/
1288 : /* ILWISOpen() */
1289 : /************************************************************************/
1290 61 : void ILWISRasterBand::ILWISOpen( string pszFileName )
1291 : {
1292 61 : ILWISDataset* dataset = (ILWISDataset*) poDS;
1293 61 : string pszDataFile;
1294 61 : pszDataFile = string(CPLResetExtension( pszFileName.c_str(), "mp#" ));
1295 :
1296 109 : fpRaw = VSIFOpenL( pszDataFile.c_str(), (dataset->eAccess == GA_Update) ? "rb+" : "rb");
1297 61 : }
1298 :
1299 : /************************************************************************/
1300 : /* ReadValueDomainProperties() */
1301 : /************************************************************************/
1302 : // Helper function for GetILWISInfo, to avoid code-duplication
1303 : // Unfortunately with side-effect (changes members psInfo and eDataType)
1304 12 : void ILWISRasterBand::ReadValueDomainProperties(string pszFileName)
1305 : {
1306 12 : string rangeString = ReadElement("BaseMap", "Range", pszFileName.c_str());
1307 24 : psInfo.vr = ValueRange(rangeString);
1308 12 : double rStep = psInfo.vr.get_rStep();
1309 12 : if ( rStep != 0 )
1310 : {
1311 10 : psInfo.bUseValueRange = true; // use ILWIS ValueRange object to convert from "raw" to "value"
1312 10 : double rMin = psInfo.vr.get_rLo();
1313 10 : double rMax = psInfo.vr.get_rHi();
1314 10 : if (rStep - (long)rStep == 0.0) // Integer values
1315 : {
1316 14 : if ( rMin >= 0 && rMax <= UCHAR_MAX)
1317 4 : eDataType = GDT_Byte;
1318 6 : else if ( rMin >= SHRT_MIN && rMax <= SHRT_MAX)
1319 0 : eDataType = GDT_Int16;
1320 6 : else if ( rMin >= 0 && rMax <= USHRT_MAX)
1321 0 : eDataType = GDT_UInt16;
1322 12 : else if ( rMin >= INT_MIN && rMax <= INT_MAX)
1323 6 : eDataType = GDT_Int32;
1324 0 : else if ( rMin >= 0 && rMax <= UINT_MAX)
1325 0 : eDataType = GDT_UInt32;
1326 : else
1327 0 : eDataType = GDT_Float64;
1328 : }
1329 : else // Floating point values
1330 : {
1331 0 : if ((rMin >= -FLT_MAX) && (rMax <= FLT_MAX) && (fabs(rStep) >= FLT_EPSILON)) // is "float" good enough?
1332 0 : eDataType = GDT_Float32;
1333 : else
1334 0 : eDataType = GDT_Float64;
1335 : }
1336 : }
1337 : else
1338 : {
1339 2 : if (psInfo.stStoreType == stFloat) // is "float" good enough?
1340 2 : eDataType = GDT_Float32;
1341 : else
1342 0 : eDataType = GDT_Float64;
1343 12 : }
1344 12 : }
1345 :
1346 : /************************************************************************/
1347 : /* GetILWISInfo() */
1348 : /************************************************************************/
1349 : // Calculates members psInfo and eDataType
1350 13 : CPLErr ILWISRasterBand::GetILWISInfo(string pszFileName)
1351 : {
1352 : // Fill the psInfo struct with defaults.
1353 : // Get the store type from the ODF
1354 13 : if (GetStoreType(pszFileName, psInfo.stStoreType) != CE_None)
1355 : {
1356 0 : return CE_Failure;
1357 : }
1358 13 : psInfo.bUseValueRange = false;
1359 13 : psInfo.stDomain = "";
1360 :
1361 : // ILWIS has several (currently 22) predefined "system-domains", that influence the data-type
1362 : // The user can also create his own domains. The possible types for these are "class", "identifier", "bool" and "value"
1363 : // The last one has Type=DomainValue
1364 : // Here we make an effort to determine the most-compact gdal-type (eDataType) that is suitable
1365 : // for the data in the current ILWIS band.
1366 : // First check against all predefined domain names (the "system-domains")
1367 : // If no match is found, read the domain ODF from disk, and get its type
1368 : // We have hardcoded the system domains here, because ILWIS may not be installed, and even if it is,
1369 : // we don't know where (thus it is useless to attempt to read a system-domain-file).
1370 :
1371 13 : string domName = ReadElement("BaseMap", "Domain", pszFileName.c_str());
1372 26 : string pszBaseName = string(CPLGetBasename( domName.c_str() ));
1373 26 : string pszPath = string(CPLGetPath( pszFileName.c_str() ));
1374 :
1375 : // Check against all "system-domains"
1376 26 : if ( EQUAL(pszBaseName.c_str(),"value") // is it a system domain with Type=DomainValue?
1377 : || EQUAL(pszBaseName.c_str(),"count")
1378 : || EQUAL(pszBaseName.c_str(),"distance")
1379 : || EQUAL(pszBaseName.c_str(),"min1to1")
1380 : || EQUAL(pszBaseName.c_str(),"nilto1")
1381 : || EQUAL(pszBaseName.c_str(),"noaa")
1382 : || EQUAL(pszBaseName.c_str(),"perc")
1383 : || EQUAL(pszBaseName.c_str(),"radar") )
1384 : {
1385 12 : ReadValueDomainProperties(pszFileName);
1386 : }
1387 1 : else if( EQUAL(pszBaseName.c_str(),"bool")
1388 : || EQUAL(pszBaseName.c_str(),"byte")
1389 : || EQUAL(pszBaseName.c_str(),"bit")
1390 : || EQUAL(pszBaseName.c_str(),"image")
1391 : || EQUAL(pszBaseName.c_str(),"colorcmp")
1392 : || EQUAL(pszBaseName.c_str(),"flowdirection")
1393 : || EQUAL(pszBaseName.c_str(),"hortonratio")
1394 : || EQUAL(pszBaseName.c_str(),"yesno") )
1395 : {
1396 0 : eDataType = GDT_Byte;
1397 0 : if( EQUAL(pszBaseName.c_str(),"image")
1398 : || EQUAL(pszBaseName.c_str(),"colorcmp"))
1399 0 : psInfo.stDomain = pszBaseName;
1400 : }
1401 1 : else if( EQUAL(pszBaseName.c_str(),"color")
1402 : || EQUAL(pszBaseName.c_str(),"none" )
1403 : || EQUAL(pszBaseName.c_str(),"coordbuf")
1404 : || EQUAL(pszBaseName.c_str(),"binary")
1405 : || EQUAL(pszBaseName.c_str(),"string") )
1406 : {
1407 : CPLError( CE_Failure, CPLE_AppDefined,
1408 0 : "Unsupported ILWIS domain type.");
1409 0 : return CE_Failure;
1410 : }
1411 : else
1412 : {
1413 : // No match found. Assume it is a self-created domain. Read its type and decide the GDAL type.
1414 1 : string pszDomainFileName = string(CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"dom" ));
1415 2 : string domType = ReadElement("Domain", "Type", pszDomainFileName.c_str());
1416 2 : if EQUAL(domType.c_str(),"domainvalue") // is it a self-created domain of type=DomainValue?
1417 : {
1418 0 : ReadValueDomainProperties(pszFileName);
1419 : }
1420 1 : else if((!EQUAL(domType.c_str(),"domainbit"))
1421 : && (!EQUAL(domType.c_str(),"domainstring"))
1422 : && (!EQUAL(domType.c_str(),"domaincolor"))
1423 : && (!EQUAL(domType.c_str(),"domainbinary"))
1424 : && (!EQUAL(domType.c_str(),"domaincoordBuf"))
1425 : && (!EQUAL(domType.c_str(),"domaincoord")))
1426 : {
1427 : // Type is "DomainClass", "DomainBool" or "DomainIdentifier".
1428 : // For now we set the GDAL storeType be the same as the ILWIS storeType
1429 : // The user will have to convert the classes manually.
1430 1 : eDataType = ILWIS2GDALType(psInfo.stStoreType);
1431 : }
1432 : else
1433 : {
1434 : CPLError( CE_Failure, CPLE_AppDefined,
1435 0 : "Unsupported ILWIS domain type.");
1436 0 : return CE_Failure;
1437 0 : }
1438 : }
1439 :
1440 13 : return CE_None;
1441 : }
1442 :
1443 : /** This driver defines a Block to be the entire raster; The method reads
1444 : each line as a block. it reads the data into pImage.
1445 :
1446 : @param nBlockXOff This must be zero for this driver
1447 : @param pImage Dump the data here
1448 :
1449 : @return A CPLErr code. This implementation returns a CE_Failure if the
1450 : block offsets are non-zero, If successful, returns CE_None. */
1451 : /************************************************************************/
1452 : /* IReadBlock() */
1453 : /************************************************************************/
1454 506 : CPLErr ILWISRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
1455 : void * pImage )
1456 : {
1457 : // pImage is empty; this function fills it with data from fpRaw
1458 : // (ILWIS data to foreign data)
1459 :
1460 : // If the x block offset is non-zero, something is wrong.
1461 : CPLAssert( nBlockXOff == 0 );
1462 :
1463 506 : int nBlockSize = nBlockXSize * nBlockYSize * nSizePerPixel;
1464 506 : if( fpRaw == NULL )
1465 : {
1466 : CPLError( CE_Failure, CPLE_OpenFailed,
1467 0 : "Failed to open ILWIS data file.");
1468 0 : return( CE_Failure );
1469 : }
1470 :
1471 : /* -------------------------------------------------------------------- */
1472 : /* Handle the case of a strip in a writable file that doesn't */
1473 : /* exist yet, but that we want to read. Just set to zeros and */
1474 : /* return. */
1475 : /* -------------------------------------------------------------------- */
1476 506 : ILWISDataset* poIDS = (ILWISDataset*) poDS;
1477 :
1478 : #ifdef notdef
1479 : if( poIDS->bNewDataset && (poIDS->eAccess == GA_Update))
1480 : {
1481 : FillWithNoData(pImage);
1482 : return CE_None;
1483 : }
1484 : #endif
1485 :
1486 506 : VSIFSeekL( fpRaw, nBlockSize*nBlockYOff, SEEK_SET );
1487 506 : void *pData = (char *)CPLMalloc(nBlockSize);
1488 506 : if (VSIFReadL( pData, 1, nBlockSize, fpRaw ) < 1)
1489 : {
1490 0 : if( poIDS->bNewDataset )
1491 : {
1492 0 : FillWithNoData(pImage);
1493 0 : return CE_None;
1494 : }
1495 : else
1496 : {
1497 0 : CPLFree( pData );
1498 : CPLError( CE_Failure, CPLE_FileIO,
1499 0 : "Read of file failed with fread error.");
1500 0 : return CE_Failure;
1501 : }
1502 : }
1503 :
1504 : // Copy the data from pData to pImage, and convert the store-type
1505 : // The data in pData has store-type = psInfo.stStoreType
1506 : // The data in pImage has store-type = eDataType
1507 : // They may not match, because we have chosen the most compact store-type,
1508 : // and for GDAL this may be different than for ILWIS.
1509 :
1510 506 : switch (psInfo.stStoreType)
1511 : {
1512 : int iCol;
1513 : case stByte:
1514 15030 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1515 : {
1516 14725 : double rV = psInfo.bUseValueRange ? psInfo.vr.rValue(((GByte *) pData)[iCol]) : ((GByte *) pData)[iCol];
1517 14725 : SetValue(pImage, iCol, rV);
1518 : }
1519 305 : break;
1520 : case stInt:
1521 0 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1522 : {
1523 0 : double rV = psInfo.bUseValueRange ? psInfo.vr.rValue(((GInt16 *) pData)[iCol]) : ((GInt16 *) pData)[iCol];
1524 0 : SetValue(pImage, iCol, rV);
1525 : }
1526 0 : break;
1527 : case stLong:
1528 0 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1529 : {
1530 0 : double rV = psInfo.bUseValueRange ? psInfo.vr.rValue(((GInt32 *) pData)[iCol]) : ((GInt32 *) pData)[iCol];
1531 0 : SetValue(pImage, iCol, rV);
1532 : }
1533 0 : break;
1534 : case stFloat:
1535 40602 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1536 40401 : ((float *) pImage)[iCol] = ((float *) pData)[iCol];
1537 201 : break;
1538 : case stReal:
1539 0 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1540 0 : ((double *) pImage)[iCol] = ((double *) pData)[iCol];
1541 : break;
1542 : default:
1543 : CPLAssert(0);
1544 : }
1545 :
1546 : // Officially we should also translate "nodata" values, but at this point
1547 : // we can't tell what's the "nodata" value of the destination (foreign) dataset
1548 :
1549 506 : CPLFree( pData );
1550 :
1551 506 : return CE_None;
1552 : }
1553 :
1554 14725 : void ILWISRasterBand::SetValue(void *pImage, int i, double rV) {
1555 14725 : switch ( eDataType ) {
1556 : case GDT_Byte:
1557 7225 : ((GByte *)pImage)[i] = (GByte) rV;
1558 7225 : break;
1559 : case GDT_UInt16:
1560 0 : ((GUInt16 *) pImage)[i] = (GUInt16) rV;
1561 0 : break;
1562 : case GDT_Int16:
1563 0 : ((GInt16 *) pImage)[i] = (GInt16) rV;
1564 0 : break;
1565 : case GDT_UInt32:
1566 0 : ((GUInt32 *) pImage)[i] = (GUInt32) rV;
1567 0 : break;
1568 : case GDT_Int32:
1569 7500 : ((GInt32 *) pImage)[i] = (GInt32) rV;
1570 7500 : break;
1571 : case GDT_Float32:
1572 0 : ((float *) pImage)[i] = (float) rV;
1573 0 : break;
1574 : case GDT_Float64:
1575 0 : ((double *) pImage)[i] = rV;
1576 : break;
1577 : default:
1578 : CPLAssert(0);
1579 : }
1580 14725 : }
1581 :
1582 7500 : double ILWISRasterBand::GetValue(void *pImage, int i) {
1583 7500 : double rV = 0; // Does GDAL have an official nodata value?
1584 7500 : switch ( eDataType ) {
1585 : case GDT_Byte:
1586 7500 : rV = ((GByte *)pImage)[i];
1587 7500 : break;
1588 : case GDT_UInt16:
1589 0 : rV = ((GUInt16 *) pImage)[i];
1590 0 : break;
1591 : case GDT_Int16:
1592 0 : rV = ((GInt16 *) pImage)[i];
1593 0 : break;
1594 : case GDT_UInt32:
1595 0 : rV = ((GUInt32 *) pImage)[i];
1596 0 : break;
1597 : case GDT_Int32:
1598 0 : rV = ((GInt32 *) pImage)[i];
1599 0 : break;
1600 : case GDT_Float32:
1601 0 : rV = ((float *) pImage)[i];
1602 0 : break;
1603 : case GDT_Float64:
1604 0 : rV = ((double *) pImage)[i];
1605 : break;
1606 : default:
1607 : CPLAssert(0);
1608 : }
1609 7500 : return rV;
1610 : }
1611 :
1612 0 : void ILWISRasterBand::FillWithNoData(void * pImage)
1613 : {
1614 0 : if (psInfo.stStoreType == stByte)
1615 0 : memset(pImage, 0, nBlockXSize * nBlockYSize);
1616 : else
1617 : {
1618 0 : switch (psInfo.stStoreType)
1619 : {
1620 : case stInt:
1621 0 : ((GInt16*)pImage)[0] = shUNDEF;
1622 0 : break;
1623 : case stLong:
1624 0 : ((GInt32*)pImage)[0] = iUNDEF;
1625 0 : break;
1626 : case stFloat:
1627 0 : ((float*)pImage)[0] = flUNDEF;
1628 0 : break;
1629 : case stReal:
1630 0 : ((double*)pImage)[0] = rUNDEF;
1631 : break;
1632 : default: // should there be handling for stByte?
1633 : break;
1634 : }
1635 0 : int iItemSize = GDALGetDataTypeSize(eDataType) / 8;
1636 0 : for (int i = 1; i < nBlockXSize * nBlockYSize; ++i)
1637 0 : memcpy( ((char*)pImage) + iItemSize * i, (char*)pImage + iItemSize * (i - 1), iItemSize);
1638 : }
1639 0 : }
1640 :
1641 : /************************************************************************/
1642 : /* IWriteBlock() */
1643 : /* */
1644 : /************************************************************************/
1645 :
1646 351 : CPLErr ILWISRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
1647 : void* pImage)
1648 : {
1649 : // pImage has data; this function reads this data and stores it to fpRaw
1650 : // (foreign data to ILWIS data)
1651 :
1652 : // Note that this function will not overwrite existing data in fpRaw, but
1653 : // it will "fill gaps" marked by "nodata" values
1654 :
1655 351 : ILWISDataset* dataset = (ILWISDataset*) poDS;
1656 :
1657 : CPLAssert( dataset != NULL
1658 : && nBlockXOff == 0
1659 : && nBlockYOff >= 0
1660 : && pImage != NULL );
1661 :
1662 351 : int nXSize = dataset->GetRasterXSize();
1663 351 : int nBlockSize = nBlockXSize * nBlockYSize * nSizePerPixel;
1664 351 : void *pData = CPLMalloc(nBlockSize);
1665 :
1666 351 : VSIFSeekL( fpRaw, nBlockSize * nBlockYOff, SEEK_SET );
1667 :
1668 351 : bool fDataExists = (VSIFReadL( pData, 1, nBlockSize, fpRaw ) >= 1);
1669 :
1670 : // Copy the data from pImage to pData, and convert the store-type
1671 : // The data in pData has store-type = psInfo.stStoreType
1672 : // The data in pImage has store-type = eDataType
1673 : // They may not match, because we have chosen the most compact store-type,
1674 : // and for GDAL this may be different than for ILWIS.
1675 :
1676 351 : if( fDataExists )
1677 : {
1678 : // fpRaw (thus pData) already has data
1679 : // Take care to not overwrite it
1680 : // thus only fill in gaps (nodata values)
1681 0 : switch (psInfo.stStoreType)
1682 : {
1683 : int iCol;
1684 : case stByte:
1685 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1686 0 : if ((( GByte * )pData)[iCol] == 0)
1687 : {
1688 0 : double rV = GetValue(pImage, iCol);
1689 0 : (( GByte * )pData)[iCol] = (GByte)
1690 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1691 : }
1692 0 : break;
1693 : case stInt:
1694 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1695 0 : if ((( GInt16 * )pData)[iCol] == shUNDEF)
1696 : {
1697 0 : double rV = GetValue(pImage, iCol);
1698 0 : (( GInt16 * )pData)[iCol] = (GInt16)
1699 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1700 : }
1701 0 : break;
1702 : case stLong:
1703 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1704 0 : if ((( GInt32 * )pData)[iCol] == iUNDEF)
1705 : {
1706 0 : double rV = GetValue(pImage, iCol);
1707 0 : (( GInt32 * )pData)[iCol] = (GInt32)
1708 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1709 : }
1710 0 : break;
1711 : case stFloat:
1712 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1713 0 : if ((( float * )pData)[iCol] == flUNDEF)
1714 0 : (( float * )pData)[iCol] = ((float* )pImage)[iCol];
1715 0 : break;
1716 : case stReal:
1717 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1718 0 : if ((( double * )pData)[iCol] == rUNDEF)
1719 0 : (( double * )pData)[iCol] = ((double* )pImage)[iCol];
1720 : break;
1721 : }
1722 : }
1723 : else
1724 : {
1725 : // fpRaw (thus pData) is still empty, just write the data
1726 351 : switch (psInfo.stStoreType)
1727 : {
1728 : int iCol;
1729 : case stByte:
1730 7650 : for (iCol = 0; iCol < nXSize; iCol++ )
1731 : {
1732 7500 : double rV = GetValue(pImage, iCol);
1733 7500 : (( GByte * )pData)[iCol] = (GByte)
1734 7500 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1735 : }
1736 150 : break;
1737 : case stInt:
1738 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1739 : {
1740 0 : double rV = GetValue(pImage, iCol);
1741 0 : (( GInt16 * )pData)[iCol] = (GInt16)
1742 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1743 : }
1744 0 : break;
1745 : case stLong:
1746 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1747 : {
1748 0 : double rV = GetValue(pImage, iCol);
1749 0 : ((GInt32 *)pData)[iCol] = (GInt32)
1750 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1751 : }
1752 0 : break;
1753 : case stFloat:
1754 40602 : for (iCol = 0; iCol < nXSize; iCol++ )
1755 40401 : (( float * )pData)[iCol] = ((float* )pImage)[iCol];
1756 201 : break;
1757 : case stReal:
1758 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1759 0 : (( double * )pData)[iCol] = ((double* )pImage)[iCol];
1760 : break;
1761 : }
1762 : }
1763 :
1764 : // Officially we should also translate "nodata" values, but at this point
1765 : // we can't tell what's the "nodata" value of the source (foreign) dataset
1766 :
1767 351 : VSIFSeekL( fpRaw, nBlockSize * nBlockYOff, SEEK_SET );
1768 :
1769 351 : if (VSIFWriteL( pData, 1, nBlockSize, fpRaw ) < 1)
1770 : {
1771 0 : CPLFree( pData );
1772 : CPLError( CE_Failure, CPLE_FileIO,
1773 0 : "Write of file failed with fwrite error.");
1774 0 : return CE_Failure;
1775 : }
1776 :
1777 351 : CPLFree( pData );
1778 351 : return CE_None;
1779 : }
1780 :
1781 : /************************************************************************/
1782 : /* GetNoDataValue() */
1783 : /************************************************************************/
1784 12 : double ILWISRasterBand::GetNoDataValue( int *pbSuccess )
1785 :
1786 : {
1787 12 : if( pbSuccess != NULL )
1788 12 : *pbSuccess = TRUE;
1789 :
1790 12 : if( eDataType == GDT_Float64 )
1791 0 : return rUNDEF;
1792 12 : else if( eDataType == GDT_Int32)
1793 3 : return iUNDEF;
1794 9 : else if( eDataType == GDT_Int16)
1795 0 : return shUNDEF;
1796 9 : else if( eDataType == GDT_Float32)
1797 2 : return flUNDEF;
1798 7 : else if( EQUAL(psInfo.stDomain.c_str(),"image")
1799 : || EQUAL(psInfo.stDomain.c_str(),"colorcmp"))
1800 : {
1801 0 : *pbSuccess = false;
1802 0 : return 0;
1803 : }
1804 : else
1805 7 : return 0;
1806 : }
1807 :
1808 : /************************************************************************/
1809 : /* ValueRange() */
1810 : /************************************************************************/
1811 0 : double Max(double a, double b)
1812 0 : { return (a>=b && a!=rUNDEF) ? a : b; }
1813 :
1814 24 : static double doubleConv(const char* s)
1815 : {
1816 24 : if (s == 0) return rUNDEF;
1817 : char *endptr;
1818 24 : char *begin = const_cast<char*>(s);
1819 :
1820 : // skip leading spaces; strtol will return 0 on a string with only spaces
1821 : // which is not what we want
1822 24 : while (isspace((unsigned char)*begin)) ++begin;
1823 :
1824 24 : if (strlen(begin) == 0) return rUNDEF;
1825 24 : errno = 0;
1826 24 : double r = strtod(begin, &endptr);
1827 24 : if ((0 == *endptr) && (errno==0))
1828 24 : return r;
1829 0 : while (*endptr != 0) { // check trailing spaces
1830 0 : if (*endptr != ' ')
1831 0 : return rUNDEF;
1832 0 : endptr++;
1833 : }
1834 0 : return r;
1835 : }
1836 :
1837 12 : ValueRange::ValueRange(string sRng)
1838 : {
1839 12 : char* sRange = new char[sRng.length() + 1];
1840 476 : for (unsigned int i = 0; i < sRng.length(); ++i)
1841 464 : sRange[i] = sRng[i];
1842 12 : sRange[sRng.length()] = 0;
1843 :
1844 12 : char *p1 = strchr(sRange, ':');
1845 12 : if (0 == p1)
1846 0 : return;
1847 :
1848 12 : char *p3 = strstr(sRange, ",offset=");
1849 12 : if (0 == p3)
1850 12 : p3 = strstr(sRange, ":offset=");
1851 12 : _r0 = rUNDEF;
1852 12 : if (0 != p3) {
1853 12 : _r0 = doubleConv(p3+8);
1854 12 : *p3 = 0;
1855 : }
1856 12 : char *p2 = strrchr(sRange, ':');
1857 12 : _rStep = 1;
1858 12 : if (p1 != p2) { // step
1859 12 : _rStep = doubleConv(p2+1);
1860 12 : *p2 = 0;
1861 : }
1862 :
1863 12 : p2 = strchr(sRange, ':');
1864 12 : if (p2 != 0) {
1865 12 : *p2 = 0;
1866 12 : _rLo = atof(sRange);
1867 12 : _rHi = atof(p2+1);
1868 : }
1869 : else {
1870 0 : _rLo = atof(sRange);
1871 0 : _rHi = _rLo;
1872 : }
1873 12 : init(_r0);
1874 :
1875 12 : delete [] sRange;
1876 : }
1877 :
1878 61 : ValueRange::ValueRange(double min, double max) // step = 1
1879 : {
1880 61 : _rLo = min;
1881 61 : _rHi = max;
1882 61 : _rStep = 1;
1883 61 : init();
1884 61 : }
1885 :
1886 0 : ValueRange::ValueRange(double min, double max, double step)
1887 : {
1888 0 : _rLo = min;
1889 0 : _rHi = max;
1890 0 : _rStep = step;
1891 0 : init();
1892 0 : }
1893 :
1894 71 : static ilwisStoreType stNeeded(unsigned long iNr)
1895 : {
1896 71 : if (iNr <= 256)
1897 65 : return stByte;
1898 6 : if (iNr <= SHRT_MAX)
1899 0 : return stInt;
1900 6 : return stLong;
1901 : }
1902 :
1903 61 : void ValueRange::init()
1904 : {
1905 61 : init(rUNDEF);
1906 61 : }
1907 :
1908 73 : void ValueRange::init(double rRaw0)
1909 : {
1910 : try {
1911 73 : _iDec = 0;
1912 73 : if (_rStep < 0)
1913 0 : _rStep = 0;
1914 73 : double r = _rStep;
1915 73 : if (r <= 1e-20)
1916 2 : _iDec = 3;
1917 142 : else while (r - floor(r) > 1e-20) {
1918 0 : r *= 10;
1919 0 : _iDec++;
1920 0 : if (_iDec > 10)
1921 0 : break;
1922 : }
1923 :
1924 73 : short iBeforeDec = 1;
1925 73 : double rMax = MAX(fabs(get_rLo()), fabs(get_rHi()));
1926 73 : if (rMax != 0)
1927 12 : iBeforeDec = (short)floor(log10(rMax)) + 1;
1928 73 : if (get_rLo() < 0)
1929 8 : iBeforeDec++;
1930 73 : _iWidth = iBeforeDec + _iDec;
1931 73 : if (_iDec > 0)
1932 2 : _iWidth++;
1933 73 : if (_iWidth > 12)
1934 0 : _iWidth = 12;
1935 73 : if (_rStep < 1e-06)
1936 : {
1937 2 : st = stReal;
1938 2 : _rStep = 0;
1939 : }
1940 : else {
1941 71 : r = get_rHi() - get_rLo();
1942 71 : if (r <= ULONG_MAX) {
1943 71 : r /= _rStep;
1944 71 : r += 1;
1945 : }
1946 71 : r += 1;
1947 71 : if (r > LONG_MAX)
1948 0 : st = stReal;
1949 : else {
1950 71 : st = stNeeded((unsigned long)floor(r+0.5));
1951 71 : if (st < stByte)
1952 0 : st = stByte;
1953 : }
1954 : }
1955 73 : if (rUNDEF != rRaw0)
1956 12 : _r0 = rRaw0;
1957 : else {
1958 61 : _r0 = 0;
1959 61 : if (st <= stByte)
1960 61 : _r0 = -1;
1961 : }
1962 73 : if (st > stInt)
1963 8 : iRawUndef = iUNDEF;
1964 65 : else if (st == stInt)
1965 0 : iRawUndef = shUNDEF;
1966 : else
1967 65 : iRawUndef = 0;
1968 : }
1969 0 : catch (std::exception*) {
1970 0 : st = stReal;
1971 0 : _r0 = 0;
1972 0 : _rStep = 0.0001;
1973 0 : _rHi = 1e300;
1974 0 : _rLo = -1e300;
1975 0 : iRawUndef = iUNDEF;
1976 : }
1977 73 : }
1978 :
1979 0 : string ValueRange::ToString()
1980 : {
1981 : char buffer[200];
1982 0 : if (fabs(get_rLo()) > 1.0e20 || fabs(get_rHi()) > 1.0e20)
1983 0 : sprintf(buffer, "%g:%g:%f:offset=%g", get_rLo(), get_rHi(), get_rStep(), get_rRaw0());
1984 0 : else if (get_iDec() >= 0)
1985 0 : sprintf(buffer, "%.*f:%.*f:%.*f:offset=%.0f", get_iDec(), get_rLo(), get_iDec(), get_rHi(), get_iDec(), get_rStep(), get_rRaw0());
1986 : else
1987 0 : sprintf(buffer, "%f:%f:%f:offset=%.0f", get_rLo(), get_rHi(), get_rStep(), get_rRaw0());
1988 0 : return string(buffer);
1989 : }
1990 :
1991 8300 : double ValueRange::rValue(int iRaw)
1992 : {
1993 8300 : if (iRaw == iUNDEF || iRaw == iRawUndef)
1994 0 : return rUNDEF;
1995 8300 : double rVal = iRaw + _r0;
1996 8300 : rVal *= _rStep;
1997 8300 : if (get_rLo() == get_rHi())
1998 0 : return rVal;
1999 8300 : double rEpsilon = _rStep == 0.0 ? 1e-6 : _rStep / 3.0; // avoid any rounding problems with an epsilon directly based on the
2000 : // the stepsize
2001 8300 : if ((rVal - get_rLo() < -rEpsilon) || (rVal - get_rHi() > rEpsilon))
2002 0 : return rUNDEF;
2003 8300 : return rVal;
2004 : }
2005 :
2006 0 : int ValueRange::iRaw(double rValue)
2007 : {
2008 0 : if (rValue == rUNDEF) // || !fContains(rValue))
2009 0 : return iUNDEF;
2010 0 : double rEpsilon = _rStep == 0.0 ? 1e-6 : _rStep / 3.0;
2011 0 : if (rValue - get_rLo() < -rEpsilon) // take a little rounding tolerance
2012 0 : return iUNDEF;
2013 0 : else if (rValue - get_rHi() > rEpsilon) // take a little rounding tolerance
2014 0 : return iUNDEF;
2015 0 : rValue /= _rStep;
2016 0 : double rVal = floor(rValue+0.5);
2017 0 : rVal -= _r0;
2018 : long iVal;
2019 0 : iVal = longConv(rVal);
2020 0 : return iVal;
2021 : }
2022 :
2023 :
2024 : /************************************************************************/
2025 : /* GDALRegister_ILWIS() */
2026 : /************************************************************************/
2027 :
2028 338 : void GDALRegister_ILWIS()
2029 :
2030 : {
2031 : GDALDriver *poDriver;
2032 :
2033 338 : if( GDALGetDriverByName( "ILWIS" ) == NULL )
2034 : {
2035 336 : poDriver = new GDALDriver();
2036 :
2037 336 : poDriver->SetDescription( "ILWIS" );
2038 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
2039 336 : "ILWIS Raster Map" );
2040 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "mpr/mpl" );
2041 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2042 336 : "Byte Int16 Int32 Float64" );
2043 336 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2044 :
2045 336 : poDriver->pfnOpen = ILWISDataset::Open;
2046 336 : poDriver->pfnCreate = ILWISDataset::Create;
2047 336 : poDriver->pfnCreateCopy = ILWISDataset::CreateCopy;
2048 :
2049 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
2050 : }
2051 338 : }
|
