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