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(FILE* 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 1808 : IniFile::IniFile(const string& filenam)
77 : {
78 1808 : filename = filenam;
79 1808 : Load();
80 1808 : bChanged = false; // Start tracking changes
81 1808 : }
82 :
83 1808 : IniFile::~IniFile()
84 : {
85 1808 : if (bChanged)
86 : {
87 1436 : Store();
88 1436 : bChanged = false;
89 : }
90 :
91 7268 : for (Sections::iterator iter = sections.begin(); iter != sections.end(); ++iter)
92 : {
93 5460 : (*(*iter).second).clear();
94 5460 : delete (*iter).second;
95 : }
96 :
97 1808 : sections.clear();
98 1808 : }
99 :
100 11916 : void IniFile::SetKeyValue(const string& section, const string& key, const string& value)
101 : {
102 11916 : Sections::iterator iterSect = sections.find(section);
103 11916 : if (iterSect == sections.end())
104 : {
105 : // Add a new section, with one new key/value entry
106 5460 : SectionEntries *entries = new SectionEntries;
107 5460 : (*entries)[key] = value;
108 5460 : 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 11916 : bChanged = true;
117 11916 : }
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 1808 : void IniFile::Load()
159 : {
160 : enum ParseState { FindSection, FindKey, ReadFindKey, StoreKey, None } state;
161 1808 : FILE *filIni = VSIFOpenL(filename.c_str(), "r");
162 1808 : if (filIni == NULL)
163 156 : 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 1436 : void IniFile::Store()
219 : {
220 1436 : FILE *filIni = VSIFOpenL(filename.c_str(), "w+");
221 1436 : if (filIni == NULL)
222 0 : 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 : bool WriteElement(string sSection, string sEntry,
273 1436 : string fn, string sValue)
274 : {
275 1436 : if (0 == fn.length())
276 0 : return false;
277 :
278 1436 : IniFile MyIniFile (fn);
279 :
280 1436 : MyIniFile.SetKeyValue(sSection, sEntry, sValue);
281 1436 : return true;
282 : }
283 :
284 : bool WriteElement(string sSection, string sEntry,
285 113 : 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 : bool WriteElement(string sSection, string sEntry,
297 186 : 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 58 : static string GDALType2ILWIS(GDALDataType type)
366 : {
367 58 : string sStoreType;
368 58 : sStoreType = "";
369 58 : switch( type )
370 : {
371 : case GDT_Byte:{
372 33 : sStoreType = "Byte";
373 33 : 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 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
544 : || 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 : double dLLLat = (adfGeoTransform[3]
552 31 : + nYSize * adfGeoTransform[5] );
553 31 : double dLLLong = (adfGeoTransform[0] );
554 31 : double dURLat = (adfGeoTransform[3] );
555 : 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 20 : CPLErr ILWISDataset::GetGeoTransform( double * padfTransform )
623 :
624 : {
625 20 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
626 20 : 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 10857 : GDALDataset *ILWISDataset::Open( GDALOpenInfo * poOpenInfo )
657 :
658 : {
659 : /* -------------------------------------------------------------------- */
660 : /* Does this look like an ILWIS file */
661 : /* -------------------------------------------------------------------- */
662 10857 : if( poOpenInfo->nHeaderBytes < 1 )
663 9581 : return NULL;
664 :
665 1276 : string sExt = CPLGetExtension( poOpenInfo->pszFilename );
666 2552 : if (!EQUAL(sExt.c_str(),"mpr") && !EQUAL(sExt.c_str(),"mpl"))
667 2542 : 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 10 : return( poDS );
816 : }
817 :
818 : /************************************************************************/
819 : /* FlushCache() */
820 : /************************************************************************/
821 :
822 59 : void ILWISDataset::FlushCache()
823 :
824 : {
825 59 : GDALDataset::FlushCache();
826 :
827 59 : if( bGeoDirty == TRUE )
828 : {
829 31 : WriteGeoReference();
830 31 : WriteProjection();
831 31 : bGeoDirty = FALSE;
832 : }
833 59 : }
834 :
835 : /************************************************************************/
836 : /* Create() */
837 : /* */
838 : /* Create a new ILWIS file. */
839 : /************************************************************************/
840 :
841 : GDALDataset *ILWISDataset::Create(const char* pszFilename,
842 : int nXSize, int nYSize,
843 : int nBands, GDALDataType eType,
844 47 : char** papszParmList)
845 : {
846 : ILWISDataset *poDS;
847 : int iBand;
848 :
849 : /* -------------------------------------------------------------------- */
850 : /* Verify input options. */
851 : /* -------------------------------------------------------------------- */
852 47 : if( eType != GDT_Byte && eType != GDT_Int16 && eType != GDT_Int32
853 : && eType != GDT_Float32 && eType != GDT_Float64 && eType != GDT_UInt16 && eType != GDT_UInt32)
854 : {
855 : CPLError( CE_Failure, CPLE_AppDefined,
856 : "Attempt to create ILWIS dataset with an illegal\n"
857 : "data type (%s).\n",
858 12 : GDALGetDataTypeName(eType) );
859 :
860 12 : return NULL;
861 : }
862 :
863 : /* -------------------------------------------------------------------- */
864 : /* Translate the data type. */
865 : /* Determine store type of ILWIS raster */
866 : /* -------------------------------------------------------------------- */
867 35 : string sDomain= "value.dom";
868 35 : double stepsize = 1;
869 35 : string sStoreType = GDALType2ILWIS(eType);
870 35 : if( EQUAL(sStoreType.c_str(),""))
871 0 : return NULL;
872 35 : else if( EQUAL(sStoreType.c_str(),"Real") || EQUAL(sStoreType.c_str(),"float"))
873 7 : stepsize = 0;
874 :
875 35 : string pszBaseName = string(CPLGetBasename( pszFilename ));
876 70 : string pszPath = string(CPLGetPath( pszFilename ));
877 :
878 : /* -------------------------------------------------------------------- */
879 : /* Write out object definition file for each band */
880 : /* -------------------------------------------------------------------- */
881 35 : string pszODFName;
882 35 : string pszDataBaseName;
883 35 : string pszFileName;
884 :
885 : char strsize[45];
886 35 : sprintf(strsize, "%d %d", nYSize, nXSize);
887 :
888 : //Form map/maplist name.
889 35 : if ( nBands == 1 )
890 : {
891 17 : pszODFName = string(CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpr"));
892 17 : pszDataBaseName = pszBaseName;
893 17 : pszFileName = CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpr");
894 : }
895 : else
896 : {
897 18 : pszFileName = CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpl");
898 18 : WriteElement("Ilwis", "Type", string(pszFileName), "MapList");
899 36 : WriteElement("MapList", "GeoRef", string(pszFileName), "none.grf");
900 36 : WriteElement("MapList", "Size", string(pszFileName), string(strsize));
901 36 : WriteElement("MapList", "Maps", string(pszFileName), nBands);
902 : }
903 :
904 35 : for( iBand = 0; iBand < nBands; iBand++ )
905 : {
906 66 : if ( nBands > 1 )
907 : {
908 : char pszBandName[100];
909 49 : sprintf(pszBandName, "%s_band_%d", pszBaseName.c_str(),iBand + 1 );
910 49 : pszODFName = string(pszBandName) + ".mpr";
911 98 : pszDataBaseName = string(pszBandName);
912 49 : sprintf(pszBandName, "Map%d", iBand);
913 98 : WriteElement("MapList", string(pszBandName), string(pszFileName), pszODFName);
914 49 : pszODFName = CPLFormFilename(pszPath.c_str(),pszDataBaseName.c_str(),"mpr");
915 : }
916 : /* -------------------------------------------------------------------- */
917 : /* Write data definition per band (.mpr) */
918 : /* -------------------------------------------------------------------- */
919 :
920 66 : WriteElement("Ilwis", "Type", pszODFName, "BaseMap");
921 132 : WriteElement("BaseMap", "Type", pszODFName, "Map");
922 132 : WriteElement("Map", "Type", pszODFName, "MapStore");
923 :
924 132 : WriteElement("BaseMap", "Domain", pszODFName, sDomain);
925 66 : string pszDataName = pszDataBaseName + ".mp#";
926 66 : WriteElement("MapStore", "Data", pszODFName, pszDataName);
927 66 : WriteElement("MapStore", "Structure", pszODFName, "Line");
928 : // sStoreType is used by ILWISRasterBand constructor to determine eDataType
929 132 : WriteElement("MapStore", "Type", pszODFName, sStoreType);
930 :
931 : // For now write-out a "Range" that is as broad as possible.
932 : // If later a better range is found (by inspecting metadata in the source dataset),
933 : // the "Range" will be overwritten by a better version.
934 : double adfMinMax[2];
935 66 : adfMinMax[0] = -9999999.9;
936 66 : adfMinMax[1] = 9999999.9;
937 : char strdouble[45];
938 66 : sprintf(strdouble, "%.3f:%.3f:%3f:offset=0", adfMinMax[0], adfMinMax[1],stepsize);
939 66 : string range = string(strdouble);
940 132 : WriteElement("BaseMap", "Range", pszODFName, range);
941 :
942 66 : WriteElement("Map", "GeoRef", pszODFName, "none.grf");
943 132 : WriteElement("Map", "Size", pszODFName, string(strsize));
944 :
945 : /* -------------------------------------------------------------------- */
946 : /* Try to create the data file. */
947 : /* -------------------------------------------------------------------- */
948 132 : pszDataName = CPLResetExtension(pszODFName.c_str(), "mp#" );
949 :
950 66 : FILE *fp = VSIFOpenL( pszDataName.c_str(), "wb" );
951 :
952 66 : if( fp == NULL )
953 : {
954 : CPLError( CE_Failure, CPLE_OpenFailed,
955 0 : "Unable to create file %s.\n", pszDataName.c_str() );
956 35 : return NULL;
957 : }
958 66 : VSIFCloseL( fp );
959 : }
960 35 : poDS = new ILWISDataset();
961 35 : poDS->nRasterXSize = nXSize;
962 35 : poDS->nRasterYSize = nYSize;
963 35 : poDS->nBands = nBands;
964 35 : poDS->eAccess = GA_Update;
965 35 : poDS->bNewDataset = TRUE;
966 35 : poDS->SetDescription(pszFilename);
967 35 : poDS->osFileName = pszFileName;
968 35 : poDS->pszIlwFileName = string(pszFileName);
969 35 : if ( nBands == 1 )
970 17 : poDS->pszFileType = "Map";
971 : else
972 18 : poDS->pszFileType = "MapList";
973 :
974 : /* -------------------------------------------------------------------- */
975 : /* Create band information objects. */
976 : /* -------------------------------------------------------------------- */
977 :
978 101 : for( iBand = 1; iBand <= poDS->nBands; iBand++ )
979 : {
980 66 : poDS->SetBand( iBand, new ILWISRasterBand( poDS, iBand ) );
981 : }
982 :
983 35 : return poDS;
984 : //return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
985 : }
986 :
987 : /************************************************************************/
988 : /* CreateCopy() */
989 : /************************************************************************/
990 :
991 : GDALDataset *
992 : ILWISDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
993 : int bStrict, char ** papszOptions,
994 18 : GDALProgressFunc pfnProgress, void * pProgressData )
995 :
996 : {
997 :
998 : ILWISDataset *poDS;
999 18 : GDALDataType eType = GDT_Byte;
1000 : int iBand;
1001 : (void) bStrict;
1002 :
1003 :
1004 18 : int nXSize = poSrcDS->GetRasterXSize();
1005 18 : int nYSize = poSrcDS->GetRasterYSize();
1006 18 : int nBands = poSrcDS->GetRasterCount();
1007 :
1008 18 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1009 0 : return NULL;
1010 :
1011 : /* -------------------------------------------------------------------- */
1012 : /* Create the basic dataset. */
1013 : /* -------------------------------------------------------------------- */
1014 45 : for( iBand = 0; iBand < nBands; iBand++ )
1015 : {
1016 27 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1017 27 : eType = GDALDataTypeUnion( eType, poBand->GetRasterDataType() );
1018 : }
1019 :
1020 : poDS = (ILWISDataset *) Create( pszFilename,
1021 : poSrcDS->GetRasterXSize(),
1022 : poSrcDS->GetRasterYSize(),
1023 : nBands,
1024 18 : eType, papszOptions );
1025 :
1026 18 : if( poDS == NULL )
1027 4 : return NULL;
1028 14 : string pszBaseName = string(CPLGetBasename( pszFilename ));
1029 28 : string pszPath = string(CPLGetPath( pszFilename ));
1030 :
1031 : /* -------------------------------------------------------------------- */
1032 : /* Copy and geo-transform and projection information. */
1033 : /* -------------------------------------------------------------------- */
1034 : double adfGeoTransform[6];
1035 28 : string georef = "none.grf";
1036 :
1037 : //check wheather we should create georeference file.
1038 : //source dataset must be north up
1039 28 : if( poSrcDS->GetGeoTransform( adfGeoTransform ) == CE_None
1040 : && (adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
1041 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
1042 : || adfGeoTransform[4] != 0.0 || fabs(adfGeoTransform[5]) != 1.0))
1043 : {
1044 13 : poDS->SetGeoTransform( adfGeoTransform );
1045 13 : if (adfGeoTransform[2] == 0.0 && adfGeoTransform[4] == 0.0)
1046 13 : georef = pszBaseName + ".grf";
1047 : }
1048 :
1049 14 : const char *pszProj = poSrcDS->GetProjectionRef();
1050 14 : if( pszProj != NULL && strlen(pszProj) > 0 )
1051 13 : poDS->SetProjection( pszProj );
1052 :
1053 : /* -------------------------------------------------------------------- */
1054 : /* Create the output raster files for each band */
1055 : /* -------------------------------------------------------------------- */
1056 :
1057 14 : for( iBand = 0; iBand < nBands; iBand++ )
1058 : {
1059 : FILE *fpData;
1060 : GByte *pData;
1061 :
1062 23 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1063 23 : ILWISRasterBand *desBand = (ILWISRasterBand *) poDS->GetRasterBand( iBand+1 );
1064 :
1065 : /* -------------------------------------------------------------------- */
1066 : /* Translate the data type. */
1067 : /* -------------------------------------------------------------------- */
1068 23 : int nLineSize = nXSize * GDALGetDataTypeSize(eType) / 8;
1069 23 : pData = (GByte *) CPLMalloc( nLineSize );
1070 :
1071 : //Determine the nodata value
1072 : int bHasNoDataValue;
1073 23 : double dNoDataValue = poBand->GetNoDataValue(&bHasNoDataValue);
1074 :
1075 : //Determine store type of ILWIS raster
1076 23 : string sStoreType = GDALType2ILWIS( eType );
1077 23 : double stepsize = 1;
1078 23 : if( EQUAL(sStoreType.c_str(),""))
1079 14 : return NULL;
1080 23 : else if( EQUAL(sStoreType.c_str(),"Real") || EQUAL(sStoreType.c_str(),"float"))
1081 2 : stepsize = 0;
1082 :
1083 : //Form the image file name, create the object definition file.
1084 23 : string pszODFName;
1085 23 : string pszDataBaseName;
1086 23 : if (nBands == 1)
1087 : {
1088 9 : pszODFName = string(CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"mpr"));
1089 9 : pszDataBaseName = pszBaseName;
1090 : }
1091 : else
1092 : {
1093 : char pszName[100];
1094 14 : sprintf(pszName, "%s_band_%d", pszBaseName.c_str(),iBand + 1 );
1095 14 : pszODFName = string(CPLFormFilename(pszPath.c_str(),pszName,"mpr"));
1096 28 : pszDataBaseName = string(pszName);
1097 : }
1098 : /* -------------------------------------------------------------------- */
1099 : /* Write data definition file for each band (.mpr) */
1100 : /* -------------------------------------------------------------------- */
1101 :
1102 : double adfMinMax[2];
1103 : int bGotMin, bGotMax;
1104 :
1105 23 : adfMinMax[0] = poBand->GetMinimum( &bGotMin );
1106 23 : adfMinMax[1] = poBand->GetMaximum( &bGotMax );
1107 23 : if( ! (bGotMin && bGotMax) )
1108 23 : GDALComputeRasterMinMax((GDALRasterBandH)poBand, FALSE, adfMinMax);
1109 23 : if ((!CPLIsNan(adfMinMax[0])) && CPLIsFinite(adfMinMax[0]) && (!CPLIsNan(adfMinMax[1])) && CPLIsFinite(adfMinMax[1]))
1110 : {
1111 : // only write a range if we got a correct one from the source dataset (otherwise ILWIS can't show the map properly)
1112 : char strdouble[45];
1113 23 : sprintf(strdouble, "%.3f:%.3f:%3f:offset=0", adfMinMax[0], adfMinMax[1],stepsize);
1114 23 : string range = string(strdouble);
1115 46 : WriteElement("BaseMap", "Range", pszODFName, range);
1116 : }
1117 23 : WriteElement("Map", "GeoRef", pszODFName, georef);
1118 :
1119 : /* -------------------------------------------------------------------- */
1120 : /* Loop over image, copy the image data. */
1121 : /* -------------------------------------------------------------------- */
1122 23 : CPLErr eErr = CE_None;
1123 :
1124 : //For file name for raw data, and create binary files.
1125 23 : string pszDataFileName = CPLResetExtension(pszODFName.c_str(), "mp#" );
1126 :
1127 23 : fpData = desBand->fpRaw;
1128 23 : if( fpData == NULL )
1129 : {
1130 : CPLError( CE_Failure, CPLE_OpenFailed,
1131 : "Attempt to create file `%s' failed.\n",
1132 0 : pszFilename );
1133 0 : return NULL;
1134 : }
1135 :
1136 273 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ )
1137 : {
1138 : eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
1139 : pData, nXSize, 1, eType,
1140 250 : 0, 0 );
1141 :
1142 250 : if( eErr == CE_None )
1143 : {
1144 250 : if (bHasNoDataValue)
1145 : {
1146 : // pData may have entries with value = dNoDataValue
1147 : // ILWIS uses a fixed value for nodata, depending on the data-type
1148 : // Therefore translate the NoDataValue from each band to ILWIS
1149 0 : for (int iCol = 0; iCol < nXSize; iCol++ )
1150 : {
1151 0 : if( EQUAL(sStoreType.c_str(),"Byte"))
1152 : {
1153 0 : if ( ((GByte * )pData)[iCol] == dNoDataValue )
1154 0 : (( GByte * )pData)[iCol] = 0;
1155 : }
1156 0 : else if( EQUAL(sStoreType.c_str(),"Int"))
1157 : {
1158 0 : if ( ((GInt16 * )pData)[iCol] == dNoDataValue )
1159 0 : (( GInt16 * )pData)[iCol] = shUNDEF;
1160 : }
1161 0 : else if( EQUAL(sStoreType.c_str(),"Long"))
1162 : {
1163 0 : if ( ((GInt32 * )pData)[iCol] == dNoDataValue )
1164 0 : (( GInt32 * )pData)[iCol] = iUNDEF;
1165 : }
1166 0 : else if( EQUAL(sStoreType.c_str(),"float"))
1167 : {
1168 0 : if ((((float * )pData)[iCol] == dNoDataValue ) || (CPLIsNan((( float * )pData)[iCol])))
1169 0 : (( float * )pData)[iCol] = flUNDEF;
1170 : }
1171 0 : else if( EQUAL(sStoreType.c_str(),"Real"))
1172 : {
1173 0 : if ((((double * )pData)[iCol] == dNoDataValue ) || (CPLIsNan((( double * )pData)[iCol])))
1174 0 : (( double * )pData)[iCol] = rUNDEF;
1175 : }
1176 : }
1177 : }
1178 250 : int iSize = VSIFWriteL( pData, 1, nLineSize, desBand->fpRaw );
1179 250 : if ( iSize < 1 )
1180 : {
1181 0 : CPLFree( pData );
1182 : //CPLFree( pData32 );
1183 : CPLError( CE_Failure, CPLE_FileIO,
1184 0 : "Write of file failed with fwrite error.");
1185 0 : return NULL;
1186 : }
1187 : }
1188 250 : if( !pfnProgress(iLine / (nYSize * nBands), NULL, pProgressData ) )
1189 0 : return NULL;
1190 : }
1191 23 : VSIFFlushL( fpData );
1192 23 : CPLFree( pData );
1193 : }
1194 :
1195 14 : poDS->FlushCache();
1196 :
1197 14 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1198 : {
1199 : CPLError( CE_Failure, CPLE_UserInterrupt,
1200 0 : "User terminated" );
1201 0 : delete poDS;
1202 :
1203 : GDALDriver *poILWISDriver =
1204 0 : (GDALDriver *) GDALGetDriverByName( "ILWIS" );
1205 0 : poILWISDriver->Delete( pszFilename );
1206 0 : return NULL;
1207 : }
1208 :
1209 14 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1210 :
1211 14 : return poDS;
1212 : }
1213 :
1214 : /************************************************************************/
1215 : /* ILWISRasterBand() */
1216 : /************************************************************************/
1217 :
1218 79 : ILWISRasterBand::ILWISRasterBand( ILWISDataset *poDS, int nBand )
1219 :
1220 : {
1221 79 : string sBandName;
1222 79 : if ( EQUAL(poDS->pszFileType.c_str(),"Map"))
1223 24 : sBandName = string(poDS->osFileName);
1224 : else //map list
1225 : {
1226 : //Form the band name
1227 : char cBandName[45];
1228 55 : sprintf( cBandName, "Map%d", nBand-1);
1229 55 : sBandName = ReadElement("MapList", string(cBandName), string(poDS->osFileName));
1230 55 : string sInputPath = string(CPLGetPath( poDS->osFileName));
1231 110 : string sBandPath = string(CPLGetPath( sBandName.c_str()));
1232 110 : string sBandBaseName = string(CPLGetBasename( sBandName.c_str()));
1233 110 : if ( 0==sBandPath.length() )
1234 55 : sBandName = string(CPLFormFilename(sInputPath.c_str(),sBandBaseName.c_str(),"mpr" ));
1235 : else
1236 0 : sBandName = string(CPLFormFilename(sBandPath.c_str(),sBandBaseName.c_str(),"mpr" ));
1237 : }
1238 :
1239 79 : if (poDS->bNewDataset)
1240 : {
1241 : // Called from Create():
1242 : // eDataType is defaulted to GDT_Byte by GDALRasterBand::GDALRasterBand
1243 : // Here we set it to match the value of sStoreType (that was set in ILWISDataset::Create)
1244 : // Unfortunately we can't take advantage of the ILWIS "ValueRange" object that would use
1245 : // the most compact storeType possible, without going through all values.
1246 66 : GetStoreType(sBandName, psInfo.stStoreType);
1247 66 : eDataType = ILWIS2GDALType(psInfo.stStoreType);
1248 : }
1249 : else // Called from Open(), thus convert ILWIS type from ODF to eDataType
1250 13 : GetILWISInfo(sBandName);
1251 79 : this->poDS = poDS;
1252 79 : this->nBand = nBand;
1253 79 : nBlockXSize = poDS->GetRasterXSize();
1254 79 : nBlockYSize = 1;
1255 79 : switch (psInfo.stStoreType)
1256 : {
1257 : case stByte:
1258 46 : nSizePerPixel = GDALGetDataTypeSize(GDT_Byte) / 8;
1259 46 : break;
1260 : case stInt:
1261 10 : nSizePerPixel = GDALGetDataTypeSize(GDT_Int16) / 8;
1262 10 : break;
1263 : case stLong:
1264 10 : nSizePerPixel = GDALGetDataTypeSize(GDT_Int32) / 8;
1265 10 : break;
1266 : case stFloat:
1267 8 : nSizePerPixel = GDALGetDataTypeSize(GDT_Float32) / 8;
1268 8 : break;
1269 : case stReal:
1270 5 : nSizePerPixel = GDALGetDataTypeSize(GDT_Float64) / 8;
1271 : break;
1272 : }
1273 79 : ILWISOpen(sBandName);
1274 79 : }
1275 :
1276 : /************************************************************************/
1277 : /* ~ILWISRasterBand() */
1278 : /************************************************************************/
1279 :
1280 79 : ILWISRasterBand::~ILWISRasterBand()
1281 :
1282 : {
1283 79 : if( fpRaw != NULL )
1284 : {
1285 79 : VSIFCloseL( fpRaw );
1286 79 : fpRaw = NULL;
1287 : }
1288 79 : }
1289 :
1290 :
1291 : /************************************************************************/
1292 : /* ILWISOpen() */
1293 : /************************************************************************/
1294 79 : void ILWISRasterBand::ILWISOpen( string pszFileName )
1295 : {
1296 79 : ILWISDataset* dataset = (ILWISDataset*) poDS;
1297 79 : string pszDataFile;
1298 79 : pszDataFile = string(CPLResetExtension( pszFileName.c_str(), "mp#" ));
1299 :
1300 145 : fpRaw = VSIFOpenL( pszDataFile.c_str(), (dataset->eAccess == GA_Update) ? "rb+" : "rb");
1301 79 : }
1302 :
1303 : /************************************************************************/
1304 : /* ReadValueDomainProperties() */
1305 : /************************************************************************/
1306 : // Helper function for GetILWISInfo, to avoid code-duplication
1307 : // Unfortunately with side-effect (changes members psInfo and eDataType)
1308 12 : void ILWISRasterBand::ReadValueDomainProperties(string pszFileName)
1309 : {
1310 12 : string rangeString = ReadElement("BaseMap", "Range", pszFileName.c_str());
1311 24 : psInfo.vr = ValueRange(rangeString);
1312 12 : double rStep = psInfo.vr.get_rStep();
1313 12 : if ( rStep != 0 )
1314 : {
1315 10 : psInfo.bUseValueRange = true; // use ILWIS ValueRange object to convert from "raw" to "value"
1316 10 : double rMin = psInfo.vr.get_rLo();
1317 10 : double rMax = psInfo.vr.get_rHi();
1318 10 : if (rStep - (long)rStep == 0.0) // Integer values
1319 : {
1320 14 : if ( rMin >= 0 && rMax <= UCHAR_MAX)
1321 4 : eDataType = GDT_Byte;
1322 6 : else if ( rMin >= SHRT_MIN && rMax <= SHRT_MAX)
1323 0 : eDataType = GDT_Int16;
1324 6 : else if ( rMin >= 0 && rMax <= USHRT_MAX)
1325 0 : eDataType = GDT_UInt16;
1326 12 : else if ( rMin >= INT_MIN && rMax <= INT_MAX)
1327 6 : eDataType = GDT_Int32;
1328 0 : else if ( rMin >= 0 && rMax <= UINT_MAX)
1329 0 : eDataType = GDT_UInt32;
1330 : else
1331 0 : eDataType = GDT_Float64;
1332 : }
1333 : else // Floating point values
1334 : {
1335 0 : if ((rMin >= -FLT_MAX) && (rMax <= FLT_MAX) && (fabs(rStep) >= FLT_EPSILON)) // is "float" good enough?
1336 0 : eDataType = GDT_Float32;
1337 : else
1338 0 : eDataType = GDT_Float64;
1339 : }
1340 : }
1341 : else
1342 : {
1343 2 : if (psInfo.stStoreType == stFloat) // is "float" good enough?
1344 2 : eDataType = GDT_Float32;
1345 : else
1346 0 : eDataType = GDT_Float64;
1347 12 : }
1348 12 : }
1349 :
1350 : /************************************************************************/
1351 : /* GetILWISInfo() */
1352 : /************************************************************************/
1353 : // Calculates members psInfo and eDataType
1354 13 : CPLErr ILWISRasterBand::GetILWISInfo(string pszFileName)
1355 : {
1356 : // Fill the psInfo struct with defaults.
1357 : // Get the store type from the ODF
1358 13 : if (GetStoreType(pszFileName, psInfo.stStoreType) != CE_None)
1359 : {
1360 0 : return CE_Failure;
1361 : }
1362 13 : psInfo.bUseValueRange = false;
1363 13 : psInfo.stDomain = "";
1364 :
1365 : // ILWIS has several (currently 22) predefined "system-domains", that influence the data-type
1366 : // The user can also create his own domains. The possible types for these are "class", "identifier", "bool" and "value"
1367 : // The last one has Type=DomainValue
1368 : // Here we make an effort to determine the most-compact gdal-type (eDataType) that is suitable
1369 : // for the data in the current ILWIS band.
1370 : // First check against all predefined domain names (the "system-domains")
1371 : // If no match is found, read the domain ODF from disk, and get its type
1372 : // We have hardcoded the system domains here, because ILWIS may not be installed, and even if it is,
1373 : // we don't know where (thus it is useless to attempt to read a system-domain-file).
1374 :
1375 13 : string domName = ReadElement("BaseMap", "Domain", pszFileName.c_str());
1376 26 : string pszBaseName = string(CPLGetBasename( domName.c_str() ));
1377 26 : string pszPath = string(CPLGetPath( pszFileName.c_str() ));
1378 :
1379 : // Check against all "system-domains"
1380 26 : if ( EQUAL(pszBaseName.c_str(),"value") // is it a system domain with Type=DomainValue?
1381 : || EQUAL(pszBaseName.c_str(),"count")
1382 : || EQUAL(pszBaseName.c_str(),"distance")
1383 : || EQUAL(pszBaseName.c_str(),"min1to1")
1384 : || EQUAL(pszBaseName.c_str(),"nilto1")
1385 : || EQUAL(pszBaseName.c_str(),"noaa")
1386 : || EQUAL(pszBaseName.c_str(),"perc")
1387 : || EQUAL(pszBaseName.c_str(),"radar") )
1388 : {
1389 12 : ReadValueDomainProperties(pszFileName);
1390 : }
1391 1 : else if( EQUAL(pszBaseName.c_str(),"bool")
1392 : || EQUAL(pszBaseName.c_str(),"byte")
1393 : || EQUAL(pszBaseName.c_str(),"bit")
1394 : || EQUAL(pszBaseName.c_str(),"image")
1395 : || EQUAL(pszBaseName.c_str(),"colorcmp")
1396 : || EQUAL(pszBaseName.c_str(),"flowdirection")
1397 : || EQUAL(pszBaseName.c_str(),"hortonratio")
1398 : || EQUAL(pszBaseName.c_str(),"yesno") )
1399 : {
1400 0 : eDataType = GDT_Byte;
1401 0 : if( EQUAL(pszBaseName.c_str(),"image")
1402 : || EQUAL(pszBaseName.c_str(),"colorcmp"))
1403 0 : psInfo.stDomain = pszBaseName;
1404 : }
1405 1 : else if( EQUAL(pszBaseName.c_str(),"color")
1406 : || EQUAL(pszBaseName.c_str(),"none" )
1407 : || EQUAL(pszBaseName.c_str(),"coordbuf")
1408 : || EQUAL(pszBaseName.c_str(),"binary")
1409 : || EQUAL(pszBaseName.c_str(),"string") )
1410 : {
1411 : CPLError( CE_Failure, CPLE_AppDefined,
1412 0 : "Unsupported ILWIS domain type.");
1413 13 : return CE_Failure;
1414 : }
1415 : else
1416 : {
1417 : // No match found. Assume it is a self-created domain. Read its type and decide the GDAL type.
1418 1 : string pszDomainFileName = string(CPLFormFilename(pszPath.c_str(),pszBaseName.c_str(),"dom" ));
1419 2 : string domType = ReadElement("Domain", "Type", pszDomainFileName.c_str());
1420 2 : if EQUAL(domType.c_str(),"domainvalue") // is it a self-created domain of type=DomainValue?
1421 : {
1422 0 : ReadValueDomainProperties(pszFileName);
1423 : }
1424 1 : else if((!EQUAL(domType.c_str(),"domainbit"))
1425 : && (!EQUAL(domType.c_str(),"domainstring"))
1426 : && (!EQUAL(domType.c_str(),"domaincolor"))
1427 : && (!EQUAL(domType.c_str(),"domainbinary"))
1428 : && (!EQUAL(domType.c_str(),"domaincoordBuf"))
1429 : && (!EQUAL(domType.c_str(),"domaincoord")))
1430 : {
1431 : // Type is "DomainClass", "DomainBool" or "DomainIdentifier".
1432 : // For now we set the GDAL storeType be the same as the ILWIS storeType
1433 : // The user will have to convert the classes manually.
1434 1 : eDataType = ILWIS2GDALType(psInfo.stStoreType);
1435 : }
1436 : else
1437 : {
1438 : CPLError( CE_Failure, CPLE_AppDefined,
1439 0 : "Unsupported ILWIS domain type.");
1440 0 : return CE_Failure;
1441 0 : }
1442 : }
1443 :
1444 13 : return CE_None;
1445 : }
1446 :
1447 : /** This driver defines a Block to be the entire raster; The method reads
1448 : each line as a block. it reads the data into pImage.
1449 :
1450 : @param nBlockXOff This must be zero for this driver
1451 : @param pImage Dump the data here
1452 :
1453 : @return A CPLErr code. This implementation returns a CE_Failure if the
1454 : block offsets are non-zero, If successful, returns CE_None. */
1455 : /************************************************************************/
1456 : /* IReadBlock() */
1457 : /************************************************************************/
1458 : CPLErr ILWISRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
1459 506 : void * pImage )
1460 : {
1461 : // pImage is empty; this function fills it with data from fpRaw
1462 : // (ILWIS data to foreign data)
1463 :
1464 : // If the x block offset is non-zero, something is wrong.
1465 506 : CPLAssert( nBlockXOff == 0 );
1466 :
1467 506 : int nBlockSize = nBlockXSize * nBlockYSize * nSizePerPixel;
1468 506 : if( fpRaw == NULL )
1469 : {
1470 : CPLError( CE_Failure, CPLE_OpenFailed,
1471 0 : "Failed to open ILWIS data file.");
1472 0 : return( CE_Failure );
1473 : }
1474 :
1475 : /* -------------------------------------------------------------------- */
1476 : /* Handle the case of a strip in a writable file that doesn't */
1477 : /* exist yet, but that we want to read. Just set to zeros and */
1478 : /* return. */
1479 : /* -------------------------------------------------------------------- */
1480 506 : ILWISDataset* poIDS = (ILWISDataset*) poDS;
1481 :
1482 : #ifdef notdef
1483 : if( poIDS->bNewDataset && (poIDS->eAccess == GA_Update))
1484 : {
1485 : FillWithNoData(pImage);
1486 : return CE_None;
1487 : }
1488 : #endif
1489 :
1490 506 : VSIFSeekL( fpRaw, nBlockSize*nBlockYOff, SEEK_SET );
1491 506 : void *pData = (char *)CPLMalloc(nBlockSize);
1492 506 : if (VSIFReadL( pData, 1, nBlockSize, fpRaw ) < 1)
1493 : {
1494 0 : if( poIDS->bNewDataset )
1495 : {
1496 0 : FillWithNoData(pImage);
1497 0 : return CE_None;
1498 : }
1499 : else
1500 : {
1501 0 : CPLFree( pData );
1502 : CPLError( CE_Failure, CPLE_FileIO,
1503 0 : "Read of file failed with fread error.");
1504 0 : return CE_Failure;
1505 : }
1506 : }
1507 :
1508 : // Copy the data from pData to pImage, and convert the store-type
1509 : // The data in pData has store-type = psInfo.stStoreType
1510 : // The data in pImage has store-type = eDataType
1511 : // They may not match, because we have chosen the most compact store-type,
1512 : // and for GDAL this may be different than for ILWIS.
1513 :
1514 506 : switch (psInfo.stStoreType)
1515 : {
1516 : int iCol;
1517 : case stByte:
1518 15030 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1519 : {
1520 14725 : double rV = psInfo.bUseValueRange ? psInfo.vr.rValue(((GByte *) pData)[iCol]) : ((GByte *) pData)[iCol];
1521 14725 : SetValue(pImage, iCol, rV);
1522 : }
1523 305 : break;
1524 : case stInt:
1525 0 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1526 : {
1527 0 : double rV = psInfo.bUseValueRange ? psInfo.vr.rValue(((GInt16 *) pData)[iCol]) : ((GInt16 *) pData)[iCol];
1528 0 : SetValue(pImage, iCol, rV);
1529 : }
1530 0 : break;
1531 : case stLong:
1532 0 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1533 : {
1534 0 : double rV = psInfo.bUseValueRange ? psInfo.vr.rValue(((GInt32 *) pData)[iCol]) : ((GInt32 *) pData)[iCol];
1535 0 : SetValue(pImage, iCol, rV);
1536 : }
1537 0 : break;
1538 : case stFloat:
1539 40602 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1540 40401 : ((float *) pImage)[iCol] = ((float *) pData)[iCol];
1541 201 : break;
1542 : case stReal:
1543 0 : for( iCol = 0; iCol < nBlockXSize; iCol++ )
1544 0 : ((double *) pImage)[iCol] = ((double *) pData)[iCol];
1545 0 : break;
1546 : default:
1547 0 : CPLAssert(0);
1548 : }
1549 :
1550 : // Officially we should also translate "nodata" values, but at this point
1551 : // we can't tell what's the "nodata" value of the destination (foreign) dataset
1552 :
1553 506 : CPLFree( pData );
1554 :
1555 506 : return CE_None;
1556 : }
1557 :
1558 14725 : void ILWISRasterBand::SetValue(void *pImage, int i, double rV) {
1559 14725 : switch ( eDataType ) {
1560 : case GDT_Byte:
1561 7225 : ((GByte *)pImage)[i] = (GByte) rV;
1562 7225 : break;
1563 : case GDT_UInt16:
1564 0 : ((GUInt16 *) pImage)[i] = (GUInt16) rV;
1565 0 : break;
1566 : case GDT_Int16:
1567 0 : ((GInt16 *) pImage)[i] = (GInt16) rV;
1568 0 : break;
1569 : case GDT_UInt32:
1570 0 : ((GUInt32 *) pImage)[i] = (GUInt32) rV;
1571 0 : break;
1572 : case GDT_Int32:
1573 7500 : ((GInt32 *) pImage)[i] = (GInt32) rV;
1574 7500 : break;
1575 : case GDT_Float32:
1576 0 : ((float *) pImage)[i] = (float) rV;
1577 0 : break;
1578 : case GDT_Float64:
1579 0 : ((double *) pImage)[i] = rV;
1580 0 : break;
1581 : default:
1582 0 : CPLAssert(0);
1583 : }
1584 14725 : }
1585 :
1586 7500 : double ILWISRasterBand::GetValue(void *pImage, int i) {
1587 7500 : double rV = 0; // Does GDAL have an official nodata value?
1588 7500 : switch ( eDataType ) {
1589 : case GDT_Byte:
1590 7500 : rV = ((GByte *)pImage)[i];
1591 7500 : break;
1592 : case GDT_UInt16:
1593 0 : rV = ((GUInt16 *) pImage)[i];
1594 0 : break;
1595 : case GDT_Int16:
1596 0 : rV = ((GInt16 *) pImage)[i];
1597 0 : break;
1598 : case GDT_UInt32:
1599 0 : rV = ((GUInt32 *) pImage)[i];
1600 0 : break;
1601 : case GDT_Int32:
1602 0 : rV = ((GInt32 *) pImage)[i];
1603 0 : break;
1604 : case GDT_Float32:
1605 0 : rV = ((float *) pImage)[i];
1606 0 : break;
1607 : case GDT_Float64:
1608 0 : rV = ((double *) pImage)[i];
1609 0 : break;
1610 : default:
1611 0 : CPLAssert(0);
1612 : }
1613 7500 : return rV;
1614 : }
1615 :
1616 0 : void ILWISRasterBand::FillWithNoData(void * pImage)
1617 : {
1618 0 : if (psInfo.stStoreType == stByte)
1619 0 : memset(pImage, 0, nBlockXSize * nBlockYSize);
1620 : else
1621 : {
1622 0 : switch (psInfo.stStoreType)
1623 : {
1624 : case stInt:
1625 0 : ((GInt16*)pImage)[0] = shUNDEF;
1626 0 : break;
1627 : case stLong:
1628 0 : ((GInt32*)pImage)[0] = iUNDEF;
1629 0 : break;
1630 : case stFloat:
1631 0 : ((float*)pImage)[0] = flUNDEF;
1632 0 : break;
1633 : case stReal:
1634 0 : ((double*)pImage)[0] = rUNDEF;
1635 : break;
1636 : default: // should there be handling for stByte?
1637 : break;
1638 : }
1639 0 : int iItemSize = GDALGetDataTypeSize(eDataType) / 8;
1640 0 : for (int i = 1; i < nBlockXSize * nBlockYSize; ++i)
1641 0 : memcpy( ((char*)pImage) + iItemSize * i, (char*)pImage + iItemSize * (i - 1), iItemSize);
1642 : }
1643 0 : }
1644 :
1645 : /************************************************************************/
1646 : /* IWriteBlock() */
1647 : /* */
1648 : /************************************************************************/
1649 :
1650 : CPLErr ILWISRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
1651 351 : void* pImage)
1652 : {
1653 : // pImage has data; this function reads this data and stores it to fpRaw
1654 : // (foreign data to ILWIS data)
1655 :
1656 : // Note that this function will not overwrite existing data in fpRaw, but
1657 : // it will "fill gaps" marked by "nodata" values
1658 :
1659 351 : ILWISDataset* dataset = (ILWISDataset*) poDS;
1660 :
1661 351 : CPLAssert( dataset != NULL
1662 : && nBlockXOff == 0
1663 : && nBlockYOff >= 0
1664 : && pImage != NULL );
1665 :
1666 351 : int nXSize = dataset->GetRasterXSize();
1667 351 : int nBlockSize = nBlockXSize * nBlockYSize * nSizePerPixel;
1668 351 : void *pData = CPLMalloc(nBlockSize);
1669 :
1670 351 : VSIFSeekL( fpRaw, nBlockSize * nBlockYOff, SEEK_SET );
1671 :
1672 351 : bool fDataExists = (VSIFReadL( pData, 1, nBlockSize, fpRaw ) >= 1);
1673 :
1674 : // Copy the data from pImage to pData, and convert the store-type
1675 : // The data in pData has store-type = psInfo.stStoreType
1676 : // The data in pImage has store-type = eDataType
1677 : // They may not match, because we have chosen the most compact store-type,
1678 : // and for GDAL this may be different than for ILWIS.
1679 :
1680 351 : if( fDataExists )
1681 : {
1682 : // fpRaw (thus pData) already has data
1683 : // Take care to not overwrite it
1684 : // thus only fill in gaps (nodata values)
1685 0 : switch (psInfo.stStoreType)
1686 : {
1687 : int iCol;
1688 : case stByte:
1689 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1690 0 : if ((( GByte * )pData)[iCol] == 0)
1691 : {
1692 0 : double rV = GetValue(pImage, iCol);
1693 : (( GByte * )pData)[iCol] = (GByte)
1694 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1695 : }
1696 0 : break;
1697 : case stInt:
1698 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1699 0 : if ((( GInt16 * )pData)[iCol] == shUNDEF)
1700 : {
1701 0 : double rV = GetValue(pImage, iCol);
1702 : (( GInt16 * )pData)[iCol] = (GInt16)
1703 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1704 : }
1705 0 : break;
1706 : case stLong:
1707 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1708 0 : if ((( GInt32 * )pData)[iCol] == iUNDEF)
1709 : {
1710 0 : double rV = GetValue(pImage, iCol);
1711 : (( GInt32 * )pData)[iCol] = (GInt32)
1712 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1713 : }
1714 0 : break;
1715 : case stFloat:
1716 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1717 0 : if ((( float * )pData)[iCol] == flUNDEF)
1718 0 : (( float * )pData)[iCol] = ((float* )pImage)[iCol];
1719 0 : break;
1720 : case stReal:
1721 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1722 0 : if ((( double * )pData)[iCol] == rUNDEF)
1723 0 : (( double * )pData)[iCol] = ((double* )pImage)[iCol];
1724 : break;
1725 : }
1726 : }
1727 : else
1728 : {
1729 : // fpRaw (thus pData) is still empty, just write the data
1730 351 : switch (psInfo.stStoreType)
1731 : {
1732 : int iCol;
1733 : case stByte:
1734 7650 : for (iCol = 0; iCol < nXSize; iCol++ )
1735 : {
1736 7500 : double rV = GetValue(pImage, iCol);
1737 : (( GByte * )pData)[iCol] = (GByte)
1738 7500 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1739 : }
1740 150 : break;
1741 : case stInt:
1742 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1743 : {
1744 0 : double rV = GetValue(pImage, iCol);
1745 : (( GInt16 * )pData)[iCol] = (GInt16)
1746 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1747 : }
1748 0 : break;
1749 : case stLong:
1750 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1751 : {
1752 0 : double rV = GetValue(pImage, iCol);
1753 : ((GInt32 *)pData)[iCol] = (GInt32)
1754 0 : (psInfo.bUseValueRange ? psInfo.vr.iRaw(rV) : rV);
1755 : }
1756 0 : break;
1757 : case stFloat:
1758 40602 : for (iCol = 0; iCol < nXSize; iCol++ )
1759 40401 : (( float * )pData)[iCol] = ((float* )pImage)[iCol];
1760 201 : break;
1761 : case stReal:
1762 0 : for (iCol = 0; iCol < nXSize; iCol++ )
1763 0 : (( double * )pData)[iCol] = ((double* )pImage)[iCol];
1764 : break;
1765 : }
1766 : }
1767 :
1768 : // Officially we should also translate "nodata" values, but at this point
1769 : // we can't tell what's the "nodata" value of the source (foreign) dataset
1770 :
1771 351 : VSIFSeekL( fpRaw, nBlockSize * nBlockYOff, SEEK_SET );
1772 :
1773 351 : if (VSIFWriteL( pData, 1, nBlockSize, fpRaw ) < 1)
1774 : {
1775 0 : CPLFree( pData );
1776 : CPLError( CE_Failure, CPLE_FileIO,
1777 0 : "Write of file failed with fwrite error.");
1778 0 : return CE_Failure;
1779 : }
1780 :
1781 351 : CPLFree( pData );
1782 351 : return CE_None;
1783 : }
1784 :
1785 : /************************************************************************/
1786 : /* GetNoDataValue() */
1787 : /************************************************************************/
1788 12 : double ILWISRasterBand::GetNoDataValue( int *pbSuccess )
1789 :
1790 : {
1791 12 : if( pbSuccess != NULL )
1792 12 : *pbSuccess = TRUE;
1793 :
1794 12 : if( eDataType == GDT_Float64 )
1795 0 : return rUNDEF;
1796 12 : else if( eDataType == GDT_Int32)
1797 3 : return iUNDEF;
1798 9 : else if( eDataType == GDT_Int16)
1799 0 : return shUNDEF;
1800 9 : else if( eDataType == GDT_Float32)
1801 2 : return flUNDEF;
1802 7 : else if( EQUAL(psInfo.stDomain.c_str(),"image")
1803 : || EQUAL(psInfo.stDomain.c_str(),"colorcmp"))
1804 : {
1805 0 : *pbSuccess = false;
1806 0 : return 0;
1807 : }
1808 : else
1809 7 : return 0;
1810 : }
1811 :
1812 : /************************************************************************/
1813 : /* ValueRange() */
1814 : /************************************************************************/
1815 0 : double Max(double a, double b)
1816 0 : { return (a>=b && a!=rUNDEF) ? a : b; }
1817 :
1818 24 : static double doubleConv(const char* s)
1819 : {
1820 24 : if (s == 0) return rUNDEF;
1821 : char *endptr;
1822 24 : char *begin = const_cast<char*>(s);
1823 :
1824 : // skip leading spaces; strtol will return 0 on a string with only spaces
1825 : // which is not what we want
1826 24 : while (isspace((unsigned char)*begin)) ++begin;
1827 :
1828 24 : if (strlen(begin) == 0) return rUNDEF;
1829 24 : errno = 0;
1830 24 : double r = strtod(begin, &endptr);
1831 24 : if ((0 == *endptr) && (errno==0))
1832 24 : return r;
1833 0 : while (*endptr != 0) { // check trailing spaces
1834 0 : if (*endptr != ' ')
1835 0 : return rUNDEF;
1836 0 : endptr++;
1837 : }
1838 0 : return r;
1839 : }
1840 :
1841 12 : ValueRange::ValueRange(string sRng)
1842 : {
1843 12 : char* sRange = new char[sRng.length() + 1];
1844 476 : for (unsigned int i = 0; i < sRng.length(); ++i)
1845 464 : sRange[i] = sRng[i];
1846 12 : sRange[sRng.length()] = 0;
1847 :
1848 12 : char *p1 = strchr(sRange, ':');
1849 12 : if (0 == p1)
1850 0 : return;
1851 :
1852 12 : char *p3 = strstr(sRange, ",offset=");
1853 12 : if (0 == p3)
1854 12 : p3 = strstr(sRange, ":offset=");
1855 12 : _r0 = rUNDEF;
1856 12 : if (0 != p3) {
1857 12 : _r0 = doubleConv(p3+8);
1858 12 : *p3 = 0;
1859 : }
1860 12 : char *p2 = strrchr(sRange, ':');
1861 12 : _rStep = 1;
1862 12 : if (p1 != p2) { // step
1863 12 : _rStep = doubleConv(p2+1);
1864 12 : *p2 = 0;
1865 : }
1866 :
1867 12 : p2 = strchr(sRange, ':');
1868 12 : if (p2 != 0) {
1869 12 : *p2 = 0;
1870 12 : _rLo = atof(sRange);
1871 12 : _rHi = atof(p2+1);
1872 : }
1873 : else {
1874 0 : _rLo = atof(sRange);
1875 0 : _rHi = _rLo;
1876 : }
1877 12 : init(_r0);
1878 :
1879 12 : delete [] sRange;
1880 : }
1881 :
1882 79 : ValueRange::ValueRange(double min, double max) // step = 1
1883 : {
1884 79 : _rLo = min;
1885 79 : _rHi = max;
1886 79 : _rStep = 1;
1887 79 : init();
1888 79 : }
1889 :
1890 0 : ValueRange::ValueRange(double min, double max, double step)
1891 : {
1892 0 : _rLo = min;
1893 0 : _rHi = max;
1894 0 : _rStep = step;
1895 0 : init();
1896 0 : }
1897 :
1898 89 : static ilwisStoreType stNeeded(unsigned long iNr)
1899 : {
1900 89 : if (iNr <= 256)
1901 83 : return stByte;
1902 6 : if (iNr <= SHRT_MAX)
1903 0 : return stInt;
1904 6 : return stLong;
1905 : }
1906 :
1907 79 : void ValueRange::init()
1908 : {
1909 79 : init(rUNDEF);
1910 79 : }
1911 :
1912 91 : void ValueRange::init(double rRaw0)
1913 : {
1914 : try {
1915 91 : _iDec = 0;
1916 91 : if (_rStep < 0)
1917 0 : _rStep = 0;
1918 91 : double r = _rStep;
1919 91 : if (r <= 1e-20)
1920 2 : _iDec = 3;
1921 178 : else while (r - floor(r) > 1e-20) {
1922 0 : r *= 10;
1923 0 : _iDec++;
1924 0 : if (_iDec > 10)
1925 0 : break;
1926 : }
1927 :
1928 91 : short iBeforeDec = 1;
1929 91 : double rMax = MAX(fabs(get_rLo()), fabs(get_rHi()));
1930 91 : if (rMax != 0)
1931 12 : iBeforeDec = (short)floor(log10(rMax)) + 1;
1932 91 : if (get_rLo() < 0)
1933 8 : iBeforeDec++;
1934 91 : _iWidth = iBeforeDec + _iDec;
1935 91 : if (_iDec > 0)
1936 2 : _iWidth++;
1937 91 : if (_iWidth > 12)
1938 0 : _iWidth = 12;
1939 91 : if (_rStep < 1e-06)
1940 : {
1941 2 : st = stReal;
1942 2 : _rStep = 0;
1943 : }
1944 : else {
1945 89 : r = get_rHi() - get_rLo();
1946 89 : if (r <= ULONG_MAX) {
1947 89 : r /= _rStep;
1948 89 : r += 1;
1949 : }
1950 89 : r += 1;
1951 89 : if (r > LONG_MAX)
1952 0 : st = stReal;
1953 : else {
1954 89 : st = stNeeded((unsigned long)floor(r+0.5));
1955 89 : if (st < stByte)
1956 0 : st = stByte;
1957 : }
1958 : }
1959 91 : if (rUNDEF != rRaw0)
1960 12 : _r0 = rRaw0;
1961 : else {
1962 79 : _r0 = 0;
1963 79 : if (st <= stByte)
1964 79 : _r0 = -1;
1965 : }
1966 91 : if (st > stInt)
1967 8 : iRawUndef = iUNDEF;
1968 83 : else if (st == stInt)
1969 0 : iRawUndef = shUNDEF;
1970 : else
1971 83 : iRawUndef = 0;
1972 : }
1973 0 : catch (std::exception*) {
1974 0 : st = stReal;
1975 0 : _r0 = 0;
1976 0 : _rStep = 0.0001;
1977 0 : _rHi = 1e300;
1978 0 : _rLo = -1e300;
1979 0 : iRawUndef = iUNDEF;
1980 : }
1981 91 : }
1982 :
1983 0 : string ValueRange::ToString()
1984 : {
1985 : char buffer[200];
1986 0 : if (fabs(get_rLo()) > 1.0e20 || fabs(get_rHi()) > 1.0e20)
1987 0 : sprintf(buffer, "%g:%g:%f:offset=%g", get_rLo(), get_rHi(), get_rStep(), get_rRaw0());
1988 0 : else if (get_iDec() >= 0)
1989 0 : sprintf(buffer, "%.*f:%.*f:%.*f:offset=%.0f", get_iDec(), get_rLo(), get_iDec(), get_rHi(), get_iDec(), get_rStep(), get_rRaw0());
1990 : else
1991 0 : sprintf(buffer, "%f:%f:%f:offset=%.0f", get_rLo(), get_rHi(), get_rStep(), get_rRaw0());
1992 0 : return string(buffer);
1993 : }
1994 :
1995 8300 : double ValueRange::rValue(int iRaw)
1996 : {
1997 8300 : if (iRaw == iUNDEF || iRaw == iRawUndef)
1998 0 : return rUNDEF;
1999 8300 : double rVal = iRaw + _r0;
2000 8300 : rVal *= _rStep;
2001 8300 : if (get_rLo() == get_rHi())
2002 0 : return rVal;
2003 8300 : double rEpsilon = _rStep == 0.0 ? 1e-6 : _rStep / 3.0; // avoid any rounding problems with an epsilon directly based on the
2004 : // the stepsize
2005 8300 : if ((rVal - get_rLo() < -rEpsilon) || (rVal - get_rHi() > rEpsilon))
2006 0 : return rUNDEF;
2007 8300 : return rVal;
2008 : }
2009 :
2010 0 : int ValueRange::iRaw(double rValue)
2011 : {
2012 0 : if (rValue == rUNDEF) // || !fContains(rValue))
2013 0 : return iUNDEF;
2014 0 : double rEpsilon = _rStep == 0.0 ? 1e-6 : _rStep / 3.0;
2015 0 : if (rValue - get_rLo() < -rEpsilon) // take a little rounding tolerance
2016 0 : return iUNDEF;
2017 0 : else if (rValue - get_rHi() > rEpsilon) // take a little rounding tolerance
2018 0 : return iUNDEF;
2019 0 : rValue /= _rStep;
2020 0 : double rVal = floor(rValue+0.5);
2021 0 : rVal -= _r0;
2022 : long iVal;
2023 0 : iVal = longConv(rVal);
2024 0 : return iVal;
2025 : }
2026 :
2027 :
2028 : /************************************************************************/
2029 : /* GDALRegister_ILWIS() */
2030 : /************************************************************************/
2031 :
2032 409 : void GDALRegister_ILWIS()
2033 :
2034 : {
2035 : GDALDriver *poDriver;
2036 :
2037 409 : if( GDALGetDriverByName( "ILWIS" ) == NULL )
2038 : {
2039 392 : poDriver = new GDALDriver();
2040 :
2041 392 : poDriver->SetDescription( "ILWIS" );
2042 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
2043 392 : "ILWIS Raster Map" );
2044 392 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "mpr/mpl" );
2045 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2046 392 : "Byte Int16 Int32 Float64" );
2047 392 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2048 :
2049 392 : poDriver->pfnOpen = ILWISDataset::Open;
2050 392 : poDriver->pfnCreate = ILWISDataset::Create;
2051 392 : poDriver->pfnCreateCopy = ILWISDataset::CreateCopy;
2052 :
2053 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
2054 : }
2055 409 : }
|
