1 : /******************************************************************************
2 : * $Id: hkvdataset.cpp 24590 2012-06-16 16:34:48Z rouault $
3 : *
4 : * Project: GView
5 : * Purpose: Implementation of Atlantis HKV labelled blob support
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2000, Frank Warmerdam
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "cpl_string.h"
32 : #include <ctype.h>
33 : #include "ogr_spatialref.h"
34 : #include "atlsci_spheroid.h"
35 :
36 : CPL_CVSID("$Id: hkvdataset.cpp 24590 2012-06-16 16:34:48Z rouault $");
37 :
38 : CPL_C_START
39 : void GDALRegister_HKV(void);
40 : CPL_C_END
41 :
42 : /************************************************************************/
43 : /* ==================================================================== */
44 : /* HKVRasterBand */
45 : /* ==================================================================== */
46 : /************************************************************************/
47 :
48 : class HKVDataset;
49 :
50 : class HKVRasterBand : public RawRasterBand
51 : {
52 : friend class HKVDataset;
53 :
54 : public:
55 : HKVRasterBand( HKVDataset *poDS, int nBand, VSILFILE * fpRaw,
56 : unsigned int nImgOffset, int nPixelOffset,
57 : int nLineOffset,
58 : GDALDataType eDataType, int bNativeOrder );
59 : virtual ~HKVRasterBand();
60 :
61 : virtual CPLErr SetNoDataValue( double );
62 : };
63 :
64 : /************************************************************************/
65 : /* HKV Spheroids */
66 : /************************************************************************/
67 :
68 : class HKVSpheroidList : public SpheroidList
69 : {
70 :
71 : public:
72 :
73 : HKVSpheroidList();
74 : ~HKVSpheroidList();
75 :
76 : };
77 :
78 41 : HKVSpheroidList :: HKVSpheroidList()
79 : {
80 41 : num_spheroids = 58;
81 41 : epsilonR = 0.1;
82 41 : epsilonI = 0.000001;
83 :
84 41 : spheroids[0].SetValuesByEqRadiusAndInvFlattening("airy-1830",6377563.396,299.3249646);
85 41 : spheroids[1].SetValuesByEqRadiusAndInvFlattening("modified-airy",6377340.189,299.3249646);
86 41 : spheroids[2].SetValuesByEqRadiusAndInvFlattening("australian-national",6378160,298.25);
87 41 : spheroids[3].SetValuesByEqRadiusAndInvFlattening("bessel-1841-namibia",6377483.865,299.1528128);
88 41 : spheroids[4].SetValuesByEqRadiusAndInvFlattening("bessel-1841",6377397.155,299.1528128);
89 41 : spheroids[5].SetValuesByEqRadiusAndInvFlattening("clarke-1858",6378294.0,294.297);
90 41 : spheroids[6].SetValuesByEqRadiusAndInvFlattening("clarke-1866",6378206.4,294.9786982);
91 41 : spheroids[7].SetValuesByEqRadiusAndInvFlattening("clarke-1880",6378249.145,293.465);
92 41 : spheroids[8].SetValuesByEqRadiusAndInvFlattening("everest-india-1830",6377276.345,300.8017);
93 41 : spheroids[9].SetValuesByEqRadiusAndInvFlattening("everest-sabah-sarawak",6377298.556,300.8017);
94 41 : spheroids[10].SetValuesByEqRadiusAndInvFlattening("everest-india-1956",6377301.243,300.8017);
95 41 : spheroids[11].SetValuesByEqRadiusAndInvFlattening("everest-malaysia-1969",6377295.664,300.8017);
96 41 : spheroids[12].SetValuesByEqRadiusAndInvFlattening("everest-malay-sing",6377304.063,300.8017);
97 41 : spheroids[13].SetValuesByEqRadiusAndInvFlattening("everest-pakistan",6377309.613,300.8017);
98 41 : spheroids[14].SetValuesByEqRadiusAndInvFlattening("modified-fisher-1960",6378155,298.3);
99 41 : spheroids[15].SetValuesByEqRadiusAndInvFlattening("helmert-1906",6378200,298.3);
100 41 : spheroids[16].SetValuesByEqRadiusAndInvFlattening("hough-1960",6378270,297);
101 41 : spheroids[17].SetValuesByEqRadiusAndInvFlattening("hughes",6378273.0,298.279);
102 41 : spheroids[18].SetValuesByEqRadiusAndInvFlattening("indonesian-1974",6378160,298.247);
103 41 : spheroids[19].SetValuesByEqRadiusAndInvFlattening("international-1924",6378388,297);
104 41 : spheroids[20].SetValuesByEqRadiusAndInvFlattening("iugc-67",6378160.0,298.254);
105 41 : spheroids[21].SetValuesByEqRadiusAndInvFlattening("iugc-75",6378140.0,298.25298);
106 41 : spheroids[22].SetValuesByEqRadiusAndInvFlattening("krassovsky-1940",6378245,298.3);
107 41 : spheroids[23].SetValuesByEqRadiusAndInvFlattening("kaula",6378165.0,292.308);
108 41 : spheroids[24].SetValuesByEqRadiusAndInvFlattening("grs-80",6378137,298.257222101);
109 41 : spheroids[25].SetValuesByEqRadiusAndInvFlattening("south-american-1969",6378160,298.25);
110 41 : spheroids[26].SetValuesByEqRadiusAndInvFlattening("wgs-72",6378135,298.26);
111 41 : spheroids[27].SetValuesByEqRadiusAndInvFlattening("wgs-84",6378137,298.257223563);
112 41 : spheroids[28].SetValuesByEqRadiusAndInvFlattening("ev-wgs-84",6378137.0,298.252841);
113 41 : spheroids[29].SetValuesByEqRadiusAndInvFlattening("ev-bessel",6377397.0,299.1976073);
114 :
115 41 : spheroids[30].SetValuesByEqRadiusAndInvFlattening("airy_1830",6377563.396,299.3249646);
116 41 : spheroids[31].SetValuesByEqRadiusAndInvFlattening("modified_airy",6377340.189,299.3249646);
117 41 : spheroids[32].SetValuesByEqRadiusAndInvFlattening("australian_national",6378160,298.25);
118 41 : spheroids[33].SetValuesByEqRadiusAndInvFlattening("bessel_1841_namibia",6377483.865,299.1528128);
119 41 : spheroids[34].SetValuesByEqRadiusAndInvFlattening("bessel_1841",6377397.155,299.1528128);
120 41 : spheroids[35].SetValuesByEqRadiusAndInvFlattening("clarke_1858",6378294.0,294.297);
121 41 : spheroids[36].SetValuesByEqRadiusAndInvFlattening("clarke_1866",6378206.4,294.9786982);
122 41 : spheroids[37].SetValuesByEqRadiusAndInvFlattening("clarke_1880",6378249.145,293.465);
123 41 : spheroids[38].SetValuesByEqRadiusAndInvFlattening("everest_india_1830",6377276.345,300.8017);
124 41 : spheroids[39].SetValuesByEqRadiusAndInvFlattening("everest_sabah_sarawak",6377298.556,300.8017);
125 41 : spheroids[40].SetValuesByEqRadiusAndInvFlattening("everest_india_1956",6377301.243,300.8017);
126 41 : spheroids[41].SetValuesByEqRadiusAndInvFlattening("everest_malaysia_1969",6377295.664,300.8017);
127 41 : spheroids[42].SetValuesByEqRadiusAndInvFlattening("everest_malay_sing",6377304.063,300.8017);
128 41 : spheroids[43].SetValuesByEqRadiusAndInvFlattening("everest_pakistan",6377309.613,300.8017);
129 41 : spheroids[44].SetValuesByEqRadiusAndInvFlattening("modified_fisher_1960",6378155,298.3);
130 41 : spheroids[45].SetValuesByEqRadiusAndInvFlattening("helmert_1906",6378200,298.3);
131 41 : spheroids[46].SetValuesByEqRadiusAndInvFlattening("hough_1960",6378270,297);
132 41 : spheroids[47].SetValuesByEqRadiusAndInvFlattening("indonesian_1974",6378160,298.247);
133 41 : spheroids[48].SetValuesByEqRadiusAndInvFlattening("international_1924",6378388,297);
134 41 : spheroids[49].SetValuesByEqRadiusAndInvFlattening("iugc_67",6378160.0,298.254);
135 41 : spheroids[50].SetValuesByEqRadiusAndInvFlattening("iugc_75",6378140.0,298.25298);
136 41 : spheroids[51].SetValuesByEqRadiusAndInvFlattening("krassovsky_1940",6378245,298.3);
137 41 : spheroids[52].SetValuesByEqRadiusAndInvFlattening("grs_80",6378137,298.257222101);
138 41 : spheroids[53].SetValuesByEqRadiusAndInvFlattening("south_american_1969",6378160,298.25);
139 41 : spheroids[54].SetValuesByEqRadiusAndInvFlattening("wgs_72",6378135,298.26);
140 41 : spheroids[55].SetValuesByEqRadiusAndInvFlattening("wgs_84",6378137,298.257223563);
141 41 : spheroids[56].SetValuesByEqRadiusAndInvFlattening("ev_wgs_84",6378137.0,298.252841);
142 41 : spheroids[57].SetValuesByEqRadiusAndInvFlattening("ev_bessel",6377397.0,299.1976073);
143 :
144 41 : }
145 :
146 41 : HKVSpheroidList::~HKVSpheroidList()
147 :
148 : {
149 41 : }
150 :
151 : CPLErr SaveHKVAttribFile( const char *pszFilenameIn,
152 : int nXSize, int nYSize, int nBands,
153 : GDALDataType eType, int bNoDataSet,
154 : double dfNoDataValue );
155 :
156 : /************************************************************************/
157 : /* ==================================================================== */
158 : /* HKVDataset */
159 : /* ==================================================================== */
160 : /************************************************************************/
161 :
162 : class HKVDataset : public RawDataset
163 : {
164 : friend class HKVRasterBand;
165 :
166 : char *pszPath;
167 : VSILFILE *fpBlob;
168 :
169 : int nGCPCount;
170 : GDAL_GCP *pasGCPList;
171 :
172 : void ProcessGeoref(const char *);
173 : void ProcessGeorefGCP(char **, const char *, double, double);
174 : void SetVersion( float version_number );
175 : float GetVersion();
176 : float MFF2version;
177 :
178 : CPLErr SetGCPProjection(const char *); /* for use in CreateCopy */
179 :
180 : GDALDataType eRasterType;
181 :
182 : void SetNoDataValue( double );
183 :
184 : char *pszProjection;
185 : char *pszGCPProjection;
186 : double adfGeoTransform[6];
187 :
188 : char **papszAttrib;
189 :
190 : int bGeorefChanged;
191 : char **papszGeoref;
192 :
193 : /* NOTE: The MFF2 format goes against GDAL's API in that nodata values are set
194 : * per-dataset rather than per-band. To compromise, for writing out, the
195 : * dataset's nodata value will be set to the last value set on any of the
196 : * raster bands.
197 : */
198 :
199 : int bNoDataSet;
200 : int bNoDataChanged;
201 : double dfNoDataValue;
202 :
203 : public:
204 : HKVDataset();
205 : virtual ~HKVDataset();
206 :
207 : virtual int GetGCPCount();
208 : virtual const char *GetGCPProjection();
209 : virtual const GDAL_GCP *GetGCPs();
210 :
211 : virtual const char *GetProjectionRef(void);
212 : virtual CPLErr GetGeoTransform( double * );
213 :
214 : virtual CPLErr SetGeoTransform( double * );
215 : virtual CPLErr SetProjection( const char * );
216 :
217 : static GDALDataset *Open( GDALOpenInfo * );
218 : static GDALDataset *Create( const char * pszFilename,
219 : int nXSize, int nYSize, int nBands,
220 : GDALDataType eType, char ** papszParmList );
221 : static GDALDataset *CreateCopy( const char * pszFilename,
222 : GDALDataset *poSrcDS,
223 : int bStrict, char ** papszOptions,
224 : GDALProgressFunc pfnProgress,
225 : void * pProgressData );
226 :
227 : static CPLErr Delete( const char * pszName );
228 : };
229 :
230 : /************************************************************************/
231 : /* ==================================================================== */
232 : /* HKVRasterBand */
233 : /* ==================================================================== */
234 : /************************************************************************/
235 :
236 : /************************************************************************/
237 : /* HKVRasterBand() */
238 : /************************************************************************/
239 :
240 91 : HKVRasterBand::HKVRasterBand( HKVDataset *poDS, int nBand, VSILFILE * fpRaw,
241 : unsigned int nImgOffset, int nPixelOffset,
242 : int nLineOffset,
243 : GDALDataType eDataType, int bNativeOrder )
244 : : RawRasterBand( (GDALDataset *) poDS, nBand,
245 : fpRaw, nImgOffset, nPixelOffset,
246 91 : nLineOffset, eDataType, bNativeOrder, TRUE )
247 :
248 : {
249 91 : this->poDS = poDS;
250 91 : this->nBand = nBand;
251 :
252 91 : nBlockXSize = poDS->GetRasterXSize();
253 91 : nBlockYSize = 1;
254 91 : }
255 :
256 : /************************************************************************/
257 : /* SetNoDataValue() */
258 : /************************************************************************/
259 :
260 0 : CPLErr HKVRasterBand::SetNoDataValue( double dfNewValue )
261 :
262 : {
263 0 : HKVDataset *poHKVDS = (HKVDataset *) poDS;
264 0 : this->RawRasterBand::SetNoDataValue( dfNewValue );
265 0 : poHKVDS->SetNoDataValue( dfNewValue );
266 :
267 0 : return CE_None;
268 : }
269 :
270 : /************************************************************************/
271 : /* ~HKVRasterBand() */
272 : /************************************************************************/
273 :
274 91 : HKVRasterBand::~HKVRasterBand()
275 :
276 : {
277 91 : }
278 :
279 : /************************************************************************/
280 : /* ==================================================================== */
281 : /* HKVDataset */
282 : /* ==================================================================== */
283 : /************************************************************************/
284 :
285 : /************************************************************************/
286 : /* HKVDataset() */
287 : /************************************************************************/
288 :
289 41 : HKVDataset::HKVDataset()
290 : {
291 41 : pszPath = NULL;
292 41 : papszAttrib = NULL;
293 41 : papszGeoref = NULL;
294 41 : bGeorefChanged = FALSE;
295 :
296 41 : nGCPCount = 0;
297 41 : pasGCPList = NULL;
298 41 : pszProjection = CPLStrdup("");
299 41 : pszGCPProjection = CPLStrdup("");
300 41 : adfGeoTransform[0] = 0.0;
301 41 : adfGeoTransform[1] = 1.0;
302 41 : adfGeoTransform[2] = 0.0;
303 41 : adfGeoTransform[3] = 0.0;
304 41 : adfGeoTransform[4] = 0.0;
305 41 : adfGeoTransform[5] = 1.0;
306 :
307 41 : bNoDataSet = FALSE;
308 41 : bNoDataChanged = FALSE;
309 :
310 : /* Initialize datasets to new version; change if necessary */
311 41 : MFF2version = (float) 1.1;
312 41 : }
313 :
314 : /************************************************************************/
315 : /* ~HKVDataset() */
316 : /************************************************************************/
317 :
318 41 : HKVDataset::~HKVDataset()
319 :
320 : {
321 41 : FlushCache();
322 41 : if( bGeorefChanged )
323 : {
324 : const char *pszFilename;
325 :
326 25 : pszFilename = CPLFormFilename(pszPath, "georef", NULL );
327 :
328 25 : CSLSave( papszGeoref, pszFilename );
329 : }
330 :
331 41 : if( bNoDataChanged )
332 : {
333 : SaveHKVAttribFile(pszPath,
334 : this->nRasterXSize,
335 : this->nRasterYSize,
336 : this->nBands,
337 : this->eRasterType,
338 : this->bNoDataSet,
339 0 : this->dfNoDataValue );
340 :
341 : }
342 :
343 41 : if( fpBlob != NULL )
344 41 : VSIFCloseL( fpBlob );
345 :
346 41 : if( nGCPCount > 0 )
347 : {
348 11 : GDALDeinitGCPs( nGCPCount, pasGCPList );
349 11 : CPLFree( pasGCPList );
350 : }
351 :
352 41 : CPLFree( pszProjection );
353 41 : CPLFree( pszGCPProjection );
354 41 : CPLFree( pszPath );
355 41 : CSLDestroy( papszGeoref );
356 41 : CSLDestroy( papszAttrib );
357 41 : }
358 :
359 : /************************************************************************/
360 : /* SetVersion() */
361 : /************************************************************************/
362 :
363 41 : void HKVDataset::SetVersion(float version_number)
364 :
365 : {
366 : //update stored info
367 41 : MFF2version = version_number;
368 41 : }
369 :
370 : /************************************************************************/
371 : /* GetVersion() */
372 : /************************************************************************/
373 :
374 0 : float HKVDataset::GetVersion()
375 :
376 : {
377 0 : return( MFF2version );
378 : }
379 :
380 : /************************************************************************/
381 : /* SetNoDataValue() */
382 : /************************************************************************/
383 :
384 0 : void HKVDataset::SetNoDataValue( double dfNewValue )
385 :
386 : {
387 :
388 0 : this->bNoDataSet = TRUE;
389 0 : this->bNoDataChanged = TRUE;
390 0 : this->dfNoDataValue = dfNewValue;
391 0 : }
392 :
393 : /************************************************************************/
394 : /* SaveHKVAttribFile() */
395 : /************************************************************************/
396 :
397 25 : CPLErr SaveHKVAttribFile( const char *pszFilenameIn,
398 : int nXSize, int nYSize, int nBands,
399 : GDALDataType eType, int bNoDataSet,
400 : double dfNoDataValue )
401 :
402 : {
403 :
404 : FILE *fp;
405 : const char *pszFilename;
406 :
407 25 : pszFilename = CPLFormFilename( pszFilenameIn, "attrib", NULL );
408 :
409 25 : fp = VSIFOpen( pszFilename, "wt" );
410 25 : if( fp == NULL )
411 : {
412 : CPLError( CE_Failure, CPLE_OpenFailed,
413 0 : "Couldn't create %s.\n", pszFilename );
414 0 : return CE_Failure;
415 : }
416 :
417 25 : fprintf( fp, "channel.enumeration = %d\n", nBands );
418 25 : fprintf( fp, "channel.interleave = { *pixel tile sequential }\n" );
419 25 : fprintf( fp, "extent.cols = %d\n", nXSize );
420 25 : fprintf( fp, "extent.rows = %d\n", nYSize );
421 :
422 25 : switch( eType )
423 : {
424 : case GDT_Byte:
425 : fprintf( fp, "pixel.encoding = "
426 10 : "{ *unsigned twos-complement ieee-754 }\n" );
427 10 : break;
428 :
429 : case GDT_UInt16:
430 : fprintf( fp, "pixel.encoding = "
431 3 : "{ *unsigned twos-complement ieee-754 }\n" );
432 3 : break;
433 :
434 : case GDT_CInt16:
435 : case GDT_Int16:
436 : fprintf( fp, "pixel.encoding = "
437 6 : "{ unsigned *twos-complement ieee-754 }\n" );
438 6 : break;
439 :
440 : case GDT_CFloat32:
441 : case GDT_Float32:
442 : fprintf( fp, "pixel.encoding = "
443 6 : "{ unsigned twos-complement *ieee-754 }\n" );
444 6 : break;
445 :
446 : default:
447 0 : CPLAssert( FALSE );
448 : }
449 :
450 25 : fprintf( fp, "pixel.size = %d\n", GDALGetDataTypeSize(eType) );
451 25 : if( GDALDataTypeIsComplex( eType ) )
452 6 : fprintf( fp, "pixel.field = { real *complex }\n" );
453 : else
454 19 : fprintf( fp, "pixel.field = { *real complex }\n" );
455 :
456 : #ifdef CPL_MSB
457 : fprintf( fp, "pixel.order = { lsbf *msbf }\n" );
458 : #else
459 25 : fprintf( fp, "pixel.order = { *lsbf msbf }\n" );
460 : #endif
461 :
462 25 : if ( bNoDataSet )
463 0 : fprintf( fp, "pixel.no_data = %f\n", dfNoDataValue );
464 :
465 : /* version information- only create the new style */
466 25 : fprintf( fp, "version = 1.1");
467 :
468 :
469 25 : VSIFClose( fp );
470 25 : return CE_None;
471 : }
472 :
473 :
474 : /************************************************************************/
475 : /* GetProjectionRef() */
476 : /************************************************************************/
477 :
478 20 : const char *HKVDataset::GetProjectionRef()
479 :
480 : {
481 20 : return( pszProjection );
482 : }
483 :
484 : /************************************************************************/
485 : /* GetGeoTransform() */
486 : /************************************************************************/
487 :
488 25 : CPLErr HKVDataset::GetGeoTransform( double * padfTransform )
489 :
490 : {
491 25 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
492 25 : return( CE_None );
493 : }
494 :
495 : /************************************************************************/
496 : /* SetGeoTransform() */
497 : /************************************************************************/
498 :
499 25 : CPLErr HKVDataset::SetGeoTransform( double * padfTransform )
500 :
501 : {
502 : char szValue[128];
503 :
504 : /* NOTE: Geotransform coordinates must match the current projection */
505 : /* of the dataset being changed (not the geotransform source). */
506 : /* ie. be in lat/longs for LL projected; UTM for UTM projected. */
507 : /* SET PROJECTION BEFORE SETTING GEOTRANSFORM TO AVOID SYNCHRONIZATION */
508 : /* PROBLEMS! */
509 :
510 : /* Update the geotransform itself */
511 25 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
512 :
513 : /* Clear previous gcps */
514 25 : if( nGCPCount > 0 )
515 : {
516 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
517 0 : CPLFree( pasGCPList );
518 : }
519 25 : nGCPCount = 0;
520 25 : pasGCPList = NULL;
521 :
522 : /* Return if the identity transform is set */
523 25 : if (adfGeoTransform[0] == 0.0 && adfGeoTransform[1] == 1.0
524 0 : && adfGeoTransform[2] == 0.0 && adfGeoTransform[3] == 0.0
525 0 : && adfGeoTransform[4] == 0.0 && adfGeoTransform[5] == 1.0 )
526 0 : return CE_None;
527 :
528 : /* Update georef text info for saving later, and */
529 : /* update GCPs to match geotransform. */
530 :
531 : double temp_lat, temp_long;
532 25 : OGRSpatialReference oUTM;
533 25 : OGRSpatialReference oLL;
534 25 : OGRCoordinateTransformation *poTransform = NULL;
535 25 : int bSuccess=TRUE;
536 : char *pszPtemp;
537 : char *pszGCPtemp;
538 :
539 : /* Projection parameter checking will have been done */
540 : /* in SetProjection. */
541 25 : if(( CSLFetchNameValue( papszGeoref, "projection.name" ) != NULL ) &&
542 : ( EQUAL(CSLFetchNameValue( papszGeoref, "projection.name" ),"UTM" )))
543 :
544 : {
545 : /* pass copies of projection info, not originals (pointers */
546 : /* get updated by importFromWkt) */
547 0 : pszPtemp = CPLStrdup(pszProjection);
548 0 : oUTM.importFromWkt(&pszPtemp);
549 0 : (oUTM.GetAttrNode("GEOGCS"))->exportToWkt(&pszGCPtemp);
550 0 : oLL.importFromWkt(&pszGCPtemp);
551 0 : poTransform = OGRCreateCoordinateTransformation( &oUTM, &oLL );
552 0 : if( poTransform == NULL )
553 : {
554 0 : bSuccess = FALSE;
555 0 : CPLErrorReset();
556 : }
557 : }
558 25 : else if ((( CSLFetchNameValue( papszGeoref, "projection.name" ) != NULL ) &&
559 : ( !EQUAL(CSLFetchNameValue( papszGeoref, "projection.name" ),"LL" ))) ||
560 : (CSLFetchNameValue( papszGeoref, "projection.name" ) == NULL ))
561 : {
562 15 : return CE_Failure;
563 : }
564 :
565 10 : nGCPCount = 0;
566 10 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5);
567 :
568 : /* -------------------------------------------------------------------- */
569 : /* top left */
570 : /* -------------------------------------------------------------------- */
571 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
572 10 : CPLFree( pasGCPList[nGCPCount].pszId );
573 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "top_left" );
574 :
575 10 : if (MFF2version > 1.0)
576 : {
577 10 : temp_lat = padfTransform[3];
578 10 : temp_long = padfTransform[0];
579 10 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
580 10 : pasGCPList[nGCPCount].dfGCPLine = 0.0;
581 : }
582 : else
583 : {
584 0 : temp_lat = padfTransform[3] + 0.5 * padfTransform[4] + 0.5 * padfTransform[5];
585 0 : temp_long = padfTransform[0] + 0.5 * padfTransform[1]+ 0.5 * padfTransform[2];
586 0 : pasGCPList[nGCPCount].dfGCPPixel = 0.5;
587 0 : pasGCPList[nGCPCount].dfGCPLine = 0.5;
588 : }
589 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
590 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
591 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
592 10 : nGCPCount++;
593 :
594 10 : if (poTransform != NULL)
595 : {
596 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
597 0 : bSuccess = FALSE;
598 : }
599 :
600 10 : if (bSuccess)
601 : {
602 10 : sprintf( szValue, "%.10f", temp_lat );
603 : papszGeoref = CSLSetNameValue( papszGeoref, "top_left.latitude",
604 10 : szValue );
605 :
606 10 : sprintf( szValue, "%.10f", temp_long );
607 : papszGeoref = CSLSetNameValue( papszGeoref, "top_left.longitude",
608 10 : szValue );
609 : }
610 :
611 : /* -------------------------------------------------------------------- */
612 : /* top_right */
613 : /* -------------------------------------------------------------------- */
614 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
615 10 : CPLFree( pasGCPList[nGCPCount].pszId );
616 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "top_right" );
617 :
618 10 : if (MFF2version > 1.0)
619 : {
620 10 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4];
621 10 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1];
622 10 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
623 10 : pasGCPList[nGCPCount].dfGCPLine = 0.0;
624 : }
625 : else
626 : {
627 0 : temp_lat = padfTransform[3] + (GetRasterXSize()-0.5) * padfTransform[4] + 0.5 * padfTransform[5];
628 0 : temp_long = padfTransform[0] + (GetRasterXSize()-0.5) * padfTransform[1] + 0.5 * padfTransform[2];
629 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()-0.5;
630 0 : pasGCPList[nGCPCount].dfGCPLine = 0.5;
631 : }
632 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
633 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
634 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
635 10 : nGCPCount++;
636 :
637 10 : if (poTransform != NULL)
638 : {
639 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
640 0 : bSuccess = FALSE;
641 : }
642 :
643 10 : if (bSuccess)
644 : {
645 10 : sprintf( szValue, "%.10f", temp_lat );
646 : papszGeoref = CSLSetNameValue( papszGeoref, "top_right.latitude",
647 10 : szValue );
648 :
649 10 : sprintf( szValue, "%.10f", temp_long );
650 : papszGeoref = CSLSetNameValue( papszGeoref, "top_right.longitude",
651 10 : szValue );
652 : }
653 :
654 : /* -------------------------------------------------------------------- */
655 : /* bottom_left */
656 : /* -------------------------------------------------------------------- */
657 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
658 10 : CPLFree( pasGCPList[nGCPCount].pszId );
659 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "bottom_left" );
660 :
661 10 : if (MFF2version > 1.0)
662 : {
663 10 : temp_lat = padfTransform[3] + GetRasterYSize() * padfTransform[5];
664 10 : temp_long = padfTransform[0] + GetRasterYSize() * padfTransform[2];
665 10 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
666 10 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize();
667 : }
668 : else
669 : {
670 0 : temp_lat = padfTransform[3] + 0.5 * padfTransform[4] + (GetRasterYSize()-0.5) * padfTransform[5];
671 0 : temp_long = padfTransform[0] + 0.5 * padfTransform[1] + (GetRasterYSize()-0.5) * padfTransform[2];
672 0 : pasGCPList[nGCPCount].dfGCPPixel = 0.5;
673 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()-0.5;
674 : }
675 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
676 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
677 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
678 10 : nGCPCount++;
679 :
680 10 : if (poTransform != NULL)
681 : {
682 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
683 0 : bSuccess = FALSE;
684 : }
685 :
686 10 : if (bSuccess)
687 : {
688 10 : sprintf( szValue, "%.10f", temp_lat );
689 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_left.latitude",
690 10 : szValue );
691 :
692 10 : sprintf( szValue, "%.10f", temp_long );
693 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_left.longitude",
694 10 : szValue );
695 : }
696 :
697 : /* -------------------------------------------------------------------- */
698 : /* bottom_right */
699 : /* -------------------------------------------------------------------- */
700 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
701 10 : CPLFree( pasGCPList[nGCPCount].pszId );
702 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "bottom_right" );
703 :
704 10 : if (MFF2version > 1.0)
705 : {
706 20 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] +
707 20 : GetRasterYSize() * padfTransform[5];
708 20 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] +
709 20 : GetRasterYSize() * padfTransform[2];
710 10 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
711 10 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize();
712 :
713 : }
714 : else
715 : {
716 0 : temp_lat = padfTransform[3] + (GetRasterXSize()-0.5) * padfTransform[4] +
717 0 : (GetRasterYSize()-0.5) * padfTransform[5];
718 0 : temp_long = padfTransform[0] + (GetRasterXSize()-0.5) * padfTransform[1] +
719 0 : (GetRasterYSize()-0.5) * padfTransform[2];
720 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()-0.5;
721 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()-0.5;
722 : }
723 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
724 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
725 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
726 10 : nGCPCount++;
727 :
728 10 : if (poTransform != NULL)
729 : {
730 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
731 0 : bSuccess = FALSE;
732 : }
733 :
734 10 : if (bSuccess)
735 : {
736 10 : sprintf( szValue, "%.10f", temp_lat );
737 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_right.latitude",
738 10 : szValue );
739 :
740 10 : sprintf( szValue, "%.10f", temp_long );
741 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_right.longitude",
742 10 : szValue );
743 : }
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Center */
747 : /* -------------------------------------------------------------------- */
748 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
749 10 : CPLFree( pasGCPList[nGCPCount].pszId );
750 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "centre" );
751 :
752 10 : if (MFF2version > 1.0)
753 : {
754 20 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] * 0.5 +
755 20 : GetRasterYSize() * padfTransform[5] * 0.5;
756 20 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] * 0.5 +
757 20 : GetRasterYSize() * padfTransform[2] * 0.5;
758 10 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()/2.0;
759 10 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()/2.0;
760 : }
761 : else
762 : {
763 0 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] * 0.5 +
764 0 : GetRasterYSize() * padfTransform[5] * 0.5;
765 0 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] * 0.5 +
766 0 : GetRasterYSize() * padfTransform[2] * 0.5;
767 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()/2.0;
768 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()/2.0;
769 : }
770 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
771 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
772 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
773 10 : nGCPCount++;
774 :
775 10 : if (poTransform != NULL)
776 : {
777 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
778 0 : bSuccess = FALSE;
779 : }
780 :
781 10 : if (bSuccess)
782 : {
783 10 : sprintf( szValue, "%.10f", temp_lat );
784 : papszGeoref = CSLSetNameValue( papszGeoref, "centre.latitude",
785 10 : szValue );
786 :
787 10 : sprintf( szValue, "%.10f", temp_long );
788 : papszGeoref = CSLSetNameValue( papszGeoref, "centre.longitude",
789 10 : szValue );
790 : }
791 :
792 10 : if (!bSuccess)
793 : {
794 : CPLError(CE_Warning,CPLE_AppDefined,
795 0 : "Warning- error setting header info in SetGeoTransform. Changes may not be saved properly.\n");
796 : }
797 :
798 10 : if (poTransform != NULL)
799 0 : delete poTransform;
800 :
801 :
802 :
803 10 : bGeorefChanged = TRUE;
804 :
805 10 : return( CE_None );
806 : }
807 :
808 10 : CPLErr HKVDataset::SetGCPProjection( const char *pszNewProjection )
809 : {
810 :
811 10 : CPLFree( pszGCPProjection );
812 10 : this->pszGCPProjection = CPLStrdup(pszNewProjection);
813 :
814 10 : return CE_None;
815 : }
816 :
817 : /************************************************************************/
818 : /* SetProjection() */
819 : /* */
820 : /* We provide very limited support for setting the projection. */
821 : /************************************************************************/
822 :
823 25 : CPLErr HKVDataset::SetProjection( const char * pszNewProjection )
824 :
825 : {
826 : HKVSpheroidList *hkvEllipsoids;
827 : double eq_radius, inv_flattening;
828 25 : OGRErr ogrerrorEq=OGRERR_NONE;
829 25 : OGRErr ogrerrorInvf=OGRERR_NONE;
830 25 : OGRErr ogrerrorOl=OGRERR_NONE;
831 :
832 25 : char *spheroid_name = NULL;
833 :
834 : /* This function is used to update a georef file */
835 :
836 :
837 : /* printf( "HKVDataset::SetProjection(%s)\n", pszNewProjection ); */
838 :
839 25 : if( !EQUALN(pszNewProjection,"GEOGCS",6)
840 : && !EQUALN(pszNewProjection,"PROJCS",6)
841 : && !EQUAL(pszNewProjection,"") )
842 : {
843 : CPLError( CE_Failure, CPLE_AppDefined,
844 : "Only OGC WKT Projections supported for writing to HKV.\n"
845 : "%s not supported.",
846 0 : pszNewProjection );
847 :
848 0 : return CE_Failure;
849 : }
850 25 : else if (EQUAL(pszNewProjection,""))
851 : {
852 0 : CPLFree( pszProjection );
853 0 : pszProjection = (char *) CPLStrdup(pszNewProjection);
854 :
855 0 : return CE_None;
856 : }
857 25 : CPLFree( pszProjection );
858 25 : pszProjection = (char *) CPLStrdup(pszNewProjection);
859 :
860 :
861 25 : OGRSpatialReference oSRS(pszNewProjection);
862 :
863 25 : if ((oSRS.GetAttrValue("PROJECTION") != NULL) &&
864 : (EQUAL(oSRS.GetAttrValue("PROJECTION"),SRS_PT_TRANSVERSE_MERCATOR)))
865 : {
866 : char *ol_txt;
867 0 : ol_txt=(char *) CPLMalloc(255);
868 0 : papszGeoref = CSLSetNameValue( papszGeoref, "projection.name", "utm" );
869 0 : sprintf(ol_txt,"%f",oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0,&ogrerrorOl));
870 : papszGeoref = CSLSetNameValue( papszGeoref, "projection.origin_longitude",
871 0 : ol_txt );
872 0 : CPLFree(ol_txt);
873 : }
874 25 : else if ((oSRS.GetAttrValue("PROJECTION") == NULL) && (oSRS.IsGeographic()))
875 : {
876 25 : papszGeoref = CSLSetNameValue( papszGeoref, "projection.name", "LL" );
877 : }
878 : else
879 : {
880 : CPLError( CE_Warning, CPLE_AppDefined,
881 0 : "Unrecognized projection.");
882 0 : return CE_Failure;
883 : }
884 25 : eq_radius = oSRS.GetSemiMajor(&ogrerrorEq);
885 25 : inv_flattening = oSRS.GetInvFlattening(&ogrerrorInvf);
886 50 : if ((ogrerrorEq == OGRERR_NONE) && (ogrerrorInvf == OGRERR_NONE))
887 : {
888 25 : hkvEllipsoids = new HKVSpheroidList;
889 50 : spheroid_name = hkvEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(eq_radius,inv_flattening);
890 25 : if (spheroid_name != NULL)
891 : {
892 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
893 25 : spheroid_name );
894 : }
895 25 : CPLFree(spheroid_name);
896 25 : delete hkvEllipsoids;
897 : }
898 : else
899 : {
900 : /* default to previous behaviour if spheroid not found by */
901 : /* radius and inverse flattening */
902 :
903 0 : if( strstr(pszNewProjection,"Bessel") != NULL )
904 : {
905 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
906 0 : "ev-bessel" );
907 : }
908 : else
909 : {
910 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
911 0 : "ev-wgs-84" );
912 : }
913 : }
914 25 : bGeorefChanged = TRUE;
915 25 : return CE_None;
916 : }
917 :
918 : /************************************************************************/
919 : /* GetGCPCount() */
920 : /************************************************************************/
921 :
922 0 : int HKVDataset::GetGCPCount()
923 :
924 : {
925 0 : return nGCPCount;
926 : }
927 :
928 : /************************************************************************/
929 : /* GetGCPProjection() */
930 : /************************************************************************/
931 :
932 0 : const char *HKVDataset::GetGCPProjection()
933 :
934 : {
935 0 : return pszGCPProjection;
936 : }
937 :
938 : /************************************************************************/
939 : /* GetGCP() */
940 : /************************************************************************/
941 :
942 0 : const GDAL_GCP *HKVDataset::GetGCPs()
943 :
944 : {
945 0 : return pasGCPList;
946 : }
947 :
948 : /************************************************************************/
949 : /* ProcessGeorefGCP() */
950 : /************************************************************************/
951 :
952 80 : void HKVDataset::ProcessGeorefGCP( char **papszGeoref, const char *pszBase,
953 : double dfRasterX, double dfRasterY )
954 :
955 : {
956 : char szFieldName[128];
957 : double dfLat, dfLong;
958 :
959 : /* -------------------------------------------------------------------- */
960 : /* Fetch the GCP from the string list. */
961 : /* -------------------------------------------------------------------- */
962 80 : sprintf( szFieldName, "%s.latitude", pszBase );
963 80 : if( CSLFetchNameValue(papszGeoref, szFieldName) == NULL )
964 75 : return;
965 : else
966 5 : dfLat = atof(CSLFetchNameValue(papszGeoref, szFieldName));
967 :
968 5 : sprintf( szFieldName, "%s.longitude", pszBase );
969 5 : if( CSLFetchNameValue(papszGeoref, szFieldName) == NULL )
970 0 : return;
971 : else
972 5 : dfLong = atof(CSLFetchNameValue(papszGeoref, szFieldName));
973 :
974 : /* -------------------------------------------------------------------- */
975 : /* Add the gcp to the internal list. */
976 : /* -------------------------------------------------------------------- */
977 5 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
978 :
979 5 : CPLFree( pasGCPList[nGCPCount].pszId );
980 :
981 5 : pasGCPList[nGCPCount].pszId = CPLStrdup( pszBase );
982 :
983 5 : pasGCPList[nGCPCount].dfGCPX = dfLong;
984 5 : pasGCPList[nGCPCount].dfGCPY = dfLat;
985 5 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
986 :
987 5 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
988 5 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
989 :
990 5 : nGCPCount++;
991 : }
992 :
993 : /************************************************************************/
994 : /* ProcessGeoref() */
995 : /************************************************************************/
996 :
997 16 : void HKVDataset::ProcessGeoref( const char * pszFilename )
998 :
999 : {
1000 : int i;
1001 16 : HKVSpheroidList *hkvEllipsoids = NULL;
1002 :
1003 : /* -------------------------------------------------------------------- */
1004 : /* Load the georef file, and boil white space away from around */
1005 : /* the equal sign. */
1006 : /* -------------------------------------------------------------------- */
1007 16 : CSLDestroy( papszGeoref );
1008 16 : papszGeoref = CSLLoad( pszFilename );
1009 16 : if( papszGeoref == NULL )
1010 0 : return;
1011 :
1012 16 : hkvEllipsoids = new HKVSpheroidList;
1013 :
1014 59 : for( i = 0; papszGeoref[i] != NULL; i++ )
1015 : {
1016 43 : int bAfterEqual = FALSE;
1017 : int iSrc, iDst;
1018 43 : char *pszLine = papszGeoref[i];
1019 :
1020 1033 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
1021 : {
1022 990 : if( bAfterEqual || pszLine[iSrc] != ' ' )
1023 : {
1024 990 : pszLine[iDst++] = pszLine[iSrc];
1025 : }
1026 :
1027 990 : if( iDst > 0 && pszLine[iDst-1] == '=' )
1028 43 : bAfterEqual = FALSE;
1029 : }
1030 43 : pszLine[iDst] = '\0';
1031 : }
1032 :
1033 : /* -------------------------------------------------------------------- */
1034 : /* Try to get GCPs, in lat/longs . */
1035 : /* -------------------------------------------------------------------- */
1036 16 : nGCPCount = 0;
1037 16 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5);
1038 :
1039 16 : if (MFF2version > 1.0)
1040 : {
1041 : ProcessGeorefGCP( papszGeoref, "top_left",
1042 16 : 0, 0 );
1043 : ProcessGeorefGCP( papszGeoref, "top_right",
1044 16 : GetRasterXSize(), 0 );
1045 : ProcessGeorefGCP( papszGeoref, "bottom_left",
1046 16 : 0, GetRasterYSize() );
1047 : ProcessGeorefGCP( papszGeoref, "bottom_right",
1048 16 : GetRasterXSize(), GetRasterYSize() );
1049 : ProcessGeorefGCP( papszGeoref, "centre",
1050 16 : GetRasterXSize()/2.0, GetRasterYSize()/2.0 );
1051 : }
1052 : else
1053 : {
1054 : ProcessGeorefGCP( papszGeoref, "top_left",
1055 0 : 0.5, 0.5 );
1056 : ProcessGeorefGCP( papszGeoref, "top_right",
1057 0 : GetRasterXSize()-0.5, 0.5 );
1058 : ProcessGeorefGCP( papszGeoref, "bottom_left",
1059 0 : 0.5, GetRasterYSize()-0.5 );
1060 : ProcessGeorefGCP( papszGeoref, "bottom_right",
1061 0 : GetRasterXSize()-0.5, GetRasterYSize()-0.5 );
1062 : ProcessGeorefGCP( papszGeoref, "centre",
1063 0 : GetRasterXSize()/2.0, GetRasterYSize()/2.0 );
1064 : }
1065 :
1066 16 : if (nGCPCount == 0)
1067 : {
1068 15 : CPLFree(pasGCPList);
1069 15 : pasGCPList = NULL;
1070 : }
1071 :
1072 : /* -------------------------------------------------------------------- */
1073 : /* Do we have a recognised projection? */
1074 : /* -------------------------------------------------------------------- */
1075 : const char *pszProjName, *pszOriginLong, *pszSpheroidName;
1076 : double eq_radius, inv_flattening;
1077 :
1078 : pszProjName = CSLFetchNameValue(papszGeoref,
1079 16 : "projection.name");
1080 : pszOriginLong = CSLFetchNameValue(papszGeoref,
1081 16 : "projection.origin_longitude");
1082 : pszSpheroidName = CSLFetchNameValue(papszGeoref,
1083 16 : "spheroid.name");
1084 :
1085 :
1086 16 : if ((pszSpheroidName != NULL) && (hkvEllipsoids->SpheroidInList(pszSpheroidName)))
1087 : {
1088 16 : eq_radius=hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
1089 16 : inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName);
1090 : }
1091 0 : else if (pszProjName != NULL)
1092 : {
1093 0 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
1094 0 : eq_radius=hkvEllipsoids->GetSpheroidEqRadius("wgs-84");
1095 0 : inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening("wgs-84");
1096 : }
1097 :
1098 17 : if( (pszProjName != NULL) && EQUAL(pszProjName,"utm") && (nGCPCount == 5) )
1099 : {
1100 : /*int nZone = (int)((atof(pszOriginLong)+184.5) / 6.0); */
1101 : int nZone;
1102 :
1103 1 : if (pszOriginLong == NULL)
1104 : {
1105 : /* If origin not specified, assume 0.0 */
1106 : CPLError(CE_Warning,CPLE_AppDefined,
1107 0 : "Warning- no projection origin longitude specified. Assuming 0.0.");
1108 0 : nZone = 31;
1109 : }
1110 : else
1111 1 : nZone = 31 + (int) floor(atof(pszOriginLong)/6.0);
1112 :
1113 1 : OGRSpatialReference oUTM;
1114 1 : OGRSpatialReference oLL;
1115 1 : OGRCoordinateTransformation *poTransform = NULL;
1116 : double dfUtmX[5], dfUtmY[5];
1117 : int gcp_index;
1118 :
1119 1 : int bSuccess = TRUE;
1120 :
1121 1 : if( pasGCPList[4].dfGCPY < 0 )
1122 0 : oUTM.SetUTM( nZone, 0 );
1123 : else
1124 1 : oUTM.SetUTM( nZone, 1 );
1125 :
1126 1 : if (pszOriginLong != NULL)
1127 : {
1128 1 : oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN,atof(pszOriginLong));
1129 1 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));
1130 : }
1131 :
1132 1 : if ((pszSpheroidName == NULL) || (EQUAL(pszSpheroidName,"wgs-84")) ||
1133 : (EQUAL(pszSpheroidName,"wgs_84")))
1134 : {
1135 0 : oUTM.SetWellKnownGeogCS( "WGS84" );
1136 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1137 : }
1138 : else
1139 : {
1140 1 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1141 : {
1142 : oUTM.SetGeogCS( "unknown","unknown",pszSpheroidName,
1143 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1144 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1145 1 : );
1146 : oLL.SetGeogCS( "unknown","unknown",pszSpheroidName,
1147 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1148 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1149 1 : );
1150 : }
1151 : else
1152 : {
1153 0 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
1154 0 : oUTM.SetWellKnownGeogCS( "WGS84" );
1155 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1156 : }
1157 : }
1158 :
1159 1 : poTransform = OGRCreateCoordinateTransformation( &oLL, &oUTM );
1160 1 : if( poTransform == NULL )
1161 : {
1162 0 : CPLErrorReset();
1163 0 : bSuccess = FALSE;
1164 : }
1165 :
1166 6 : for(gcp_index=0;gcp_index<5;gcp_index++)
1167 : {
1168 5 : dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
1169 5 : dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
1170 :
1171 5 : if( bSuccess && !poTransform->Transform( 1, &(dfUtmX[gcp_index]), &(dfUtmY[gcp_index]) ) )
1172 0 : bSuccess = FALSE;
1173 :
1174 : }
1175 :
1176 1 : if( bSuccess )
1177 : {
1178 1 : int transform_ok = FALSE;
1179 :
1180 : /* update GCPS to proper projection */
1181 6 : for(gcp_index=0;gcp_index<5;gcp_index++)
1182 : {
1183 5 : pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
1184 5 : pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
1185 : }
1186 :
1187 1 : CPLFree( pszGCPProjection );
1188 1 : pszGCPProjection = NULL;
1189 1 : oUTM.exportToWkt( &pszGCPProjection );
1190 :
1191 1 : transform_ok = GDALGCPsToGeoTransform(5,pasGCPList,adfGeoTransform,0);
1192 :
1193 1 : CPLFree( pszProjection );
1194 1 : pszProjection = NULL;
1195 1 : if (transform_ok == FALSE)
1196 : {
1197 : /* transform may not be sufficient in all cases (slant range projection) */
1198 0 : adfGeoTransform[0] = 0.0;
1199 0 : adfGeoTransform[1] = 1.0;
1200 0 : adfGeoTransform[2] = 0.0;
1201 0 : adfGeoTransform[3] = 0.0;
1202 0 : adfGeoTransform[4] = 0.0;
1203 0 : adfGeoTransform[5] = 1.0;
1204 0 : pszProjection = CPLStrdup("");
1205 : }
1206 : else
1207 1 : oUTM.exportToWkt( &pszProjection );
1208 :
1209 : }
1210 :
1211 1 : if( poTransform != NULL )
1212 1 : delete poTransform;
1213 : }
1214 15 : else if ((pszProjName != NULL) && (nGCPCount == 5))
1215 : {
1216 0 : OGRSpatialReference oLL;
1217 0 : int transform_ok = FALSE;
1218 :
1219 :
1220 0 : if (pszOriginLong != NULL)
1221 : {
1222 0 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));
1223 : }
1224 :
1225 0 : if ((pszSpheroidName == NULL) || (EQUAL(pszSpheroidName,"wgs-84")) ||
1226 : (EQUAL(pszSpheroidName,"wgs_84")))
1227 : {
1228 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1229 : }
1230 : else
1231 : {
1232 0 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1233 : {
1234 : oLL.SetGeogCS( "","",pszSpheroidName,
1235 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1236 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1237 0 : );
1238 : }
1239 : else
1240 : {
1241 0 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
1242 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1243 : }
1244 : }
1245 :
1246 0 : transform_ok = GDALGCPsToGeoTransform(5,pasGCPList,adfGeoTransform,0);
1247 :
1248 0 : CPLFree( pszProjection );
1249 0 : pszProjection = NULL;
1250 :
1251 0 : if (transform_ok == FALSE)
1252 : {
1253 0 : adfGeoTransform[0] = 0.0;
1254 0 : adfGeoTransform[1] = 1.0;
1255 0 : adfGeoTransform[2] = 0.0;
1256 0 : adfGeoTransform[3] = 0.0;
1257 0 : adfGeoTransform[4] = 0.0;
1258 0 : adfGeoTransform[5] = 1.0;
1259 : }
1260 : else
1261 : {
1262 0 : oLL.exportToWkt( &pszProjection );
1263 : }
1264 :
1265 0 : CPLFree( pszGCPProjection );
1266 0 : pszGCPProjection = NULL;
1267 0 : oLL.exportToWkt( &pszGCPProjection );
1268 :
1269 : }
1270 :
1271 16 : delete hkvEllipsoids;
1272 : }
1273 :
1274 : /************************************************************************/
1275 : /* Open() */
1276 : /************************************************************************/
1277 :
1278 12083 : GDALDataset *HKVDataset::Open( GDALOpenInfo * poOpenInfo )
1279 :
1280 : {
1281 12083 : int i, bNoDataSet = FALSE;
1282 12083 : double dfNoDataValue = 0.0;
1283 : char **papszAttrib;
1284 : const char *pszFilename, *pszValue;
1285 : VSIStatBuf sStat;
1286 :
1287 : /* -------------------------------------------------------------------- */
1288 : /* We assume the dataset is passed as a directory. Check for */
1289 : /* an attrib and blob file as a minimum. */
1290 : /* -------------------------------------------------------------------- */
1291 12083 : if( !poOpenInfo->bIsDirectory )
1292 12041 : return NULL;
1293 :
1294 42 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "image_data", NULL);
1295 42 : if( VSIStat(pszFilename,&sStat) != 0 )
1296 1 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "blob", NULL );
1297 42 : if( VSIStat(pszFilename,&sStat) != 0 )
1298 1 : return NULL;
1299 :
1300 41 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "attrib", NULL );
1301 41 : if( VSIStat(pszFilename,&sStat) != 0 )
1302 0 : return NULL;
1303 :
1304 : /* -------------------------------------------------------------------- */
1305 : /* Load the attrib file, and boil white space away from around */
1306 : /* the equal sign. */
1307 : /* -------------------------------------------------------------------- */
1308 41 : papszAttrib = CSLLoad( pszFilename );
1309 41 : if( papszAttrib == NULL )
1310 0 : return NULL;
1311 :
1312 410 : for( i = 0; papszAttrib[i] != NULL; i++ )
1313 : {
1314 369 : int bAfterEqual = FALSE;
1315 : int iSrc, iDst;
1316 369 : char *pszLine = papszAttrib[i];
1317 :
1318 10417 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
1319 : {
1320 10048 : if( bAfterEqual || pszLine[iSrc] != ' ' )
1321 : {
1322 8736 : pszLine[iDst++] = pszLine[iSrc];
1323 : }
1324 :
1325 10048 : if( iDst > 0 && pszLine[iDst-1] == '=' )
1326 738 : bAfterEqual = FALSE;
1327 : }
1328 369 : pszLine[iDst] = '\0';
1329 : }
1330 :
1331 : /* -------------------------------------------------------------------- */
1332 : /* Create a corresponding GDALDataset. */
1333 : /* -------------------------------------------------------------------- */
1334 : HKVDataset *poDS;
1335 :
1336 41 : poDS = new HKVDataset();
1337 :
1338 41 : poDS->pszPath = CPLStrdup( poOpenInfo->pszFilename );
1339 41 : poDS->papszAttrib = papszAttrib;
1340 :
1341 41 : poDS->eAccess = poOpenInfo->eAccess;
1342 :
1343 : /* -------------------------------------------------------------------- */
1344 : /* Set some dataset wide information. */
1345 : /* -------------------------------------------------------------------- */
1346 : int bNative, bComplex;
1347 41 : int nRawBands = 0;
1348 :
1349 82 : if( CSLFetchNameValue( papszAttrib, "extent.cols" ) == NULL
1350 : || CSLFetchNameValue( papszAttrib, "extent.rows" ) == NULL )
1351 0 : return NULL;
1352 :
1353 41 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib,"extent.cols"));
1354 41 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib,"extent.rows"));
1355 :
1356 41 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
1357 : {
1358 0 : delete poDS;
1359 0 : return NULL;
1360 : }
1361 :
1362 41 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.order");
1363 41 : if( pszValue == NULL )
1364 0 : bNative = TRUE;
1365 : else
1366 : {
1367 : #ifdef CPL_MSB
1368 : bNative = (strstr(pszValue,"*msbf") != NULL);
1369 : #else
1370 41 : bNative = (strstr(pszValue,"*lsbf") != NULL);
1371 : #endif
1372 : }
1373 :
1374 41 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.no_data");
1375 41 : if( pszValue != NULL )
1376 : {
1377 0 : bNoDataSet = TRUE;
1378 0 : dfNoDataValue = atof(pszValue);
1379 : }
1380 :
1381 41 : pszValue = CSLFetchNameValue(papszAttrib,"channel.enumeration");
1382 41 : if( pszValue != NULL )
1383 41 : nRawBands = atoi(pszValue);
1384 : else
1385 0 : nRawBands = 1;
1386 :
1387 41 : if (!GDALCheckBandCount(nRawBands, TRUE))
1388 : {
1389 0 : delete poDS;
1390 0 : return NULL;
1391 : }
1392 :
1393 41 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.field");
1394 51 : if( pszValue != NULL && strstr(pszValue,"*complex") != NULL )
1395 10 : bComplex = TRUE;
1396 : else
1397 31 : bComplex = FALSE;
1398 :
1399 : /* Get the version number, if present (if not, assume old version. */
1400 : /* Versions differ in their interpretation of corner coordinates. */
1401 :
1402 41 : if (CSLFetchNameValue( papszAttrib, "version" ) != NULL)
1403 : poDS->SetVersion((float)
1404 41 : atof(CSLFetchNameValue(papszAttrib, "version")));
1405 : else
1406 0 : poDS->SetVersion(1.0);
1407 :
1408 : /* -------------------------------------------------------------------- */
1409 : /* Figure out the datatype */
1410 : /* -------------------------------------------------------------------- */
1411 : const char * pszEncoding;
1412 41 : int nSize = 1;
1413 : int nPseudoBands;
1414 : GDALDataType eType;
1415 :
1416 41 : pszEncoding = CSLFetchNameValue(papszAttrib,"pixel.encoding");
1417 41 : if( pszEncoding == NULL )
1418 0 : pszEncoding = "{ *unsigned }";
1419 :
1420 41 : if( CSLFetchNameValue(papszAttrib,"pixel.size") != NULL )
1421 41 : nSize = atoi(CSLFetchNameValue(papszAttrib,"pixel.size"))/8;
1422 :
1423 41 : if( bComplex )
1424 10 : nPseudoBands = 2;
1425 : else
1426 31 : nPseudoBands = 1;
1427 :
1428 41 : if( nSize == 1 )
1429 16 : eType = GDT_Byte;
1430 30 : else if( nSize == 2 && strstr(pszEncoding,"*unsigned") != NULL )
1431 5 : eType = GDT_UInt16;
1432 25 : else if( nSize == 4 && bComplex )
1433 5 : eType = GDT_CInt16;
1434 15 : else if( nSize == 2 )
1435 5 : eType = GDT_Int16;
1436 10 : else if( nSize == 4 && strstr(pszEncoding,"*unsigned") != NULL )
1437 0 : eType = GDT_UInt32;
1438 10 : else if( nSize == 8 && strstr(pszEncoding,"*two") != NULL && bComplex )
1439 0 : eType = GDT_CInt32;
1440 10 : else if( nSize == 4 && strstr(pszEncoding,"*two") != NULL )
1441 0 : eType = GDT_Int32;
1442 15 : else if( nSize == 8 && bComplex )
1443 5 : eType = GDT_CFloat32;
1444 5 : else if( nSize == 4 )
1445 5 : eType = GDT_Float32;
1446 0 : else if( nSize == 16 && bComplex )
1447 0 : eType = GDT_CFloat64;
1448 0 : else if( nSize == 8 )
1449 0 : eType = GDT_Float64;
1450 : else
1451 : {
1452 : CPLError( CE_Failure, CPLE_AppDefined,
1453 : "Unsupported pixel data type in %s.\n"
1454 : "pixel.size=%d pixel.encoding=%s\n",
1455 0 : poDS->pszPath, nSize, pszEncoding );
1456 0 : delete poDS;
1457 0 : return NULL;
1458 : }
1459 :
1460 : /* -------------------------------------------------------------------- */
1461 : /* Open the blob file. */
1462 : /* -------------------------------------------------------------------- */
1463 41 : pszFilename = CPLFormFilename(poDS->pszPath, "image_data", NULL );
1464 41 : if( VSIStat(pszFilename,&sStat) != 0 )
1465 0 : pszFilename = CPLFormFilename(poDS->pszPath, "blob", NULL );
1466 41 : if( poOpenInfo->eAccess == GA_ReadOnly )
1467 : {
1468 16 : poDS->fpBlob = VSIFOpenL( pszFilename, "rb" );
1469 16 : if( poDS->fpBlob == NULL )
1470 : {
1471 : CPLError( CE_Failure, CPLE_OpenFailed,
1472 : "Unable to open file %s for read access.\n",
1473 0 : pszFilename );
1474 0 : delete poDS;
1475 0 : return NULL;
1476 : }
1477 : }
1478 : else
1479 : {
1480 25 : poDS->fpBlob = VSIFOpenL( pszFilename, "rb+" );
1481 25 : if( poDS->fpBlob == NULL )
1482 : {
1483 : CPLError( CE_Failure, CPLE_OpenFailed,
1484 : "Unable to open file %s for update access.\n",
1485 0 : pszFilename );
1486 0 : delete poDS;
1487 0 : return NULL;
1488 : }
1489 : }
1490 :
1491 : /* -------------------------------------------------------------------- */
1492 : /* Build the overview filename, as blob file = "_ovr". */
1493 : /* -------------------------------------------------------------------- */
1494 41 : char *pszOvrFilename = (char *) CPLMalloc(strlen(pszFilename)+5);
1495 :
1496 41 : sprintf( pszOvrFilename, "%s_ovr", pszFilename );
1497 :
1498 : /* -------------------------------------------------------------------- */
1499 : /* Define the bands. */
1500 : /* -------------------------------------------------------------------- */
1501 : int nPixelOffset, nLineOffset, nOffset;
1502 :
1503 41 : nPixelOffset = nRawBands * nSize;
1504 41 : nLineOffset = nPixelOffset * poDS->GetRasterXSize();
1505 41 : nOffset = 0;
1506 :
1507 132 : for( int iRawBand=0; iRawBand < nRawBands; iRawBand++ )
1508 : {
1509 : HKVRasterBand *poBand;
1510 :
1511 : poBand =
1512 : new HKVRasterBand( poDS, poDS->GetRasterCount()+1, poDS->fpBlob,
1513 : nOffset, nPixelOffset, nLineOffset,
1514 91 : eType, bNative );
1515 91 : poDS->SetBand( poDS->GetRasterCount()+1, poBand );
1516 91 : nOffset += GDALGetDataTypeSize( eType ) / 8;
1517 :
1518 91 : if( bNoDataSet )
1519 0 : poBand->SetNoDataValue( dfNoDataValue );
1520 : }
1521 :
1522 41 : poDS->eRasterType = eType;
1523 :
1524 : /* -------------------------------------------------------------------- */
1525 : /* Process the georef file if there is one. */
1526 : /* -------------------------------------------------------------------- */
1527 41 : pszFilename = CPLFormFilename(poDS->pszPath, "georef", NULL );
1528 41 : if( VSIStat(pszFilename,&sStat) == 0 )
1529 16 : poDS->ProcessGeoref(pszFilename);
1530 :
1531 : /* -------------------------------------------------------------------- */
1532 : /* Initialize any PAM information. */
1533 : /* -------------------------------------------------------------------- */
1534 41 : poDS->SetDescription( pszOvrFilename );
1535 41 : poDS->TryLoadXML();
1536 :
1537 : /* -------------------------------------------------------------------- */
1538 : /* Handle overviews. */
1539 : /* -------------------------------------------------------------------- */
1540 41 : poDS->oOvManager.Initialize( poDS, pszOvrFilename, NULL, TRUE );
1541 :
1542 41 : CPLFree( pszOvrFilename );
1543 :
1544 41 : return( poDS );
1545 : }
1546 :
1547 : /************************************************************************/
1548 : /* Create() */
1549 : /************************************************************************/
1550 :
1551 53 : GDALDataset *HKVDataset::Create( const char * pszFilenameIn,
1552 : int nXSize, int nYSize, int nBands,
1553 : GDALDataType eType,
1554 : char ** /* papszParmList */ )
1555 :
1556 : {
1557 : /* -------------------------------------------------------------------- */
1558 : /* Verify input options. */
1559 : /* -------------------------------------------------------------------- */
1560 53 : if (nBands <= 0)
1561 : {
1562 : CPLError( CE_Failure, CPLE_NotSupported,
1563 1 : "HKV driver does not support %d bands.\n", nBands);
1564 1 : return NULL;
1565 : }
1566 :
1567 52 : if( eType != GDT_Byte
1568 : && eType != GDT_UInt16 && eType != GDT_Int16
1569 : && eType != GDT_CInt16 && eType != GDT_Float32
1570 : && eType != GDT_CFloat32 )
1571 : {
1572 : CPLError( CE_Failure, CPLE_AppDefined,
1573 : "Attempt to create HKV file with currently unsupported\n"
1574 : "data type (%s).\n",
1575 15 : GDALGetDataTypeName(eType) );
1576 :
1577 15 : return NULL;
1578 : }
1579 :
1580 : /* -------------------------------------------------------------------- */
1581 : /* Establish the name of the directory we will be creating the */
1582 : /* new HKV directory in. Verify that this is a directory. */
1583 : /* -------------------------------------------------------------------- */
1584 : char *pszBaseDir;
1585 : VSIStatBuf sStat;
1586 :
1587 37 : if( strlen(CPLGetPath(pszFilenameIn)) == 0 )
1588 0 : pszBaseDir = CPLStrdup(".");
1589 : else
1590 37 : pszBaseDir = CPLStrdup(CPLGetPath(pszFilenameIn));
1591 :
1592 37 : if( CPLStat( pszBaseDir, &sStat ) != 0 || !VSI_ISDIR( sStat.st_mode ) )
1593 : {
1594 : CPLError( CE_Failure, CPLE_AppDefined,
1595 : "Attempt to create HKV dataset under %s,\n"
1596 : "but this is not a valid directory.\n",
1597 12 : pszBaseDir);
1598 12 : CPLFree( pszBaseDir );
1599 12 : return NULL;
1600 : }
1601 :
1602 25 : CPLFree( pszBaseDir );
1603 25 : pszBaseDir = NULL;
1604 :
1605 25 : if( VSIMkdir( pszFilenameIn, 0755 ) != 0 )
1606 : {
1607 : CPLError( CE_Failure, CPLE_AppDefined,
1608 : "Unable to create directory %s.\n",
1609 0 : pszFilenameIn );
1610 0 : return NULL;
1611 : }
1612 :
1613 : /* -------------------------------------------------------------------- */
1614 : /* Create the header file. */
1615 : /* -------------------------------------------------------------------- */
1616 : CPLErr CEHeaderCreated;
1617 :
1618 : CEHeaderCreated = SaveHKVAttribFile( pszFilenameIn, nXSize, nYSize,
1619 25 : nBands, eType, FALSE, 0.0 );
1620 :
1621 25 : if (CEHeaderCreated != CE_None )
1622 0 : return NULL;
1623 :
1624 : /* -------------------------------------------------------------------- */
1625 : /* Create the blob file. */
1626 : /* -------------------------------------------------------------------- */
1627 :
1628 : FILE *fp;
1629 : const char *pszFilename;
1630 :
1631 25 : pszFilename = CPLFormFilename( pszFilenameIn, "image_data", NULL );
1632 25 : fp = VSIFOpen( pszFilename, "wb" );
1633 25 : if( fp == NULL )
1634 : {
1635 : CPLError( CE_Failure, CPLE_OpenFailed,
1636 0 : "Couldn't create %s.\n", pszFilename );
1637 0 : return NULL;
1638 : }
1639 :
1640 25 : VSIFWrite( (void*)"", 1, 1, fp );
1641 25 : VSIFClose( fp );
1642 :
1643 : /* -------------------------------------------------------------------- */
1644 : /* Open the dataset normally. */
1645 : /* -------------------------------------------------------------------- */
1646 25 : return (GDALDataset *) GDALOpen( pszFilenameIn, GA_Update );
1647 : }
1648 :
1649 : /************************************************************************/
1650 : /* Delete() */
1651 : /* */
1652 : /* An HKV Blob dataset consists of a bunch of files in a */
1653 : /* directory. Try to delete all the files, then the */
1654 : /* directory. */
1655 : /************************************************************************/
1656 :
1657 0 : CPLErr HKVDataset::Delete( const char * pszName )
1658 :
1659 : {
1660 : VSIStatBuf sStat;
1661 : char **papszFiles;
1662 : int i;
1663 :
1664 0 : if( CPLStat( pszName, &sStat ) != 0
1665 : || !VSI_ISDIR(sStat.st_mode) )
1666 : {
1667 : CPLError( CE_Failure, CPLE_AppDefined,
1668 : "%s does not appear to be an HKV Dataset, as it is not\n"
1669 : "a path to a directory.",
1670 0 : pszName );
1671 0 : return CE_Failure;
1672 : }
1673 :
1674 0 : papszFiles = CPLReadDir( pszName );
1675 0 : for( i = 0; i < CSLCount(papszFiles); i++ )
1676 : {
1677 : const char *pszTarget;
1678 :
1679 0 : if( EQUAL(papszFiles[i],".") || EQUAL(papszFiles[i],"..") )
1680 0 : continue;
1681 :
1682 0 : pszTarget = CPLFormFilename(pszName, papszFiles[i], NULL );
1683 0 : if( VSIUnlink(pszTarget) != 0 )
1684 : {
1685 : CPLError( CE_Failure, CPLE_AppDefined,
1686 : "Unable to delete file %s,\n"
1687 : "HKVDataset Delete(%s) failed.\n",
1688 : pszTarget,
1689 0 : pszName );
1690 0 : CSLDestroy( papszFiles );
1691 0 : return CE_Failure;
1692 : }
1693 : }
1694 :
1695 0 : CSLDestroy( papszFiles );
1696 :
1697 0 : if( VSIRmdir( pszName ) != 0 )
1698 : {
1699 : CPLError( CE_Failure, CPLE_AppDefined,
1700 : "Unable to delete directory %s,\n"
1701 : "HKVDataset Delete() failed.\n",
1702 0 : pszName );
1703 0 : return CE_Failure;
1704 : }
1705 :
1706 0 : return CE_None;
1707 : }
1708 :
1709 : /************************************************************************/
1710 : /* CreateCopy() */
1711 : /************************************************************************/
1712 :
1713 : GDALDataset *
1714 28 : HKVDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
1715 : int bStrict, char ** papszOptions,
1716 : GDALProgressFunc pfnProgress, void * pProgressData )
1717 :
1718 : {
1719 : HKVDataset *poDS;
1720 : GDALDataType eType;
1721 : int iBand;
1722 :
1723 28 : int nBands = poSrcDS->GetRasterCount();
1724 28 : if (nBands == 0)
1725 : {
1726 : CPLError( CE_Failure, CPLE_NotSupported,
1727 1 : "HKV driver does not support source dataset with zero band.\n");
1728 1 : return NULL;
1729 : }
1730 :
1731 27 : eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1732 :
1733 27 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1734 0 : return NULL;
1735 :
1736 : /* check that other bands match type- sets type */
1737 : /* to unknown if they differ. */
1738 47 : for( iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++ )
1739 : {
1740 20 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1741 20 : eType = GDALDataTypeUnion( eType, poBand->GetRasterDataType() );
1742 : }
1743 :
1744 : poDS = (HKVDataset *) Create( pszFilename,
1745 : poSrcDS->GetRasterXSize(),
1746 : poSrcDS->GetRasterYSize(),
1747 : poSrcDS->GetRasterCount(),
1748 27 : eType, papszOptions );
1749 :
1750 : /* Check that Create worked- return Null if it didn't */
1751 27 : if (poDS == NULL)
1752 17 : return NULL;
1753 :
1754 : /* -------------------------------------------------------------------- */
1755 : /* Copy the image data. */
1756 : /* -------------------------------------------------------------------- */
1757 10 : int nXSize = poDS->GetRasterXSize();
1758 10 : int nYSize = poDS->GetRasterYSize();
1759 : int nBlockXSize, nBlockYSize, nBlockTotal, nBlocksDone;
1760 :
1761 10 : poDS->GetRasterBand(1)->GetBlockSize( &nBlockXSize, &nBlockYSize );
1762 :
1763 : nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize)
1764 : * ((nYSize + nBlockYSize - 1) / nBlockYSize)
1765 10 : * poSrcDS->GetRasterCount();
1766 :
1767 10 : nBlocksDone = 0;
1768 30 : for( iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++ )
1769 : {
1770 20 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand+1 );
1771 20 : GDALRasterBand *poDstBand = poDS->GetRasterBand( iBand+1 );
1772 : int iYOffset, iXOffset;
1773 : void *pData;
1774 : CPLErr eErr;
1775 : int pbSuccess;
1776 20 : double dfSrcNoDataValue =0.0;
1777 :
1778 : /* Get nodata value, if relevant */
1779 20 : dfSrcNoDataValue = poSrcBand->GetNoDataValue( &pbSuccess );
1780 20 : if ( pbSuccess )
1781 0 : poDS->SetNoDataValue( dfSrcNoDataValue );
1782 :
1783 : pData = CPLMalloc(nBlockXSize * nBlockYSize
1784 20 : * GDALGetDataTypeSize(eType) / 8);
1785 :
1786 220 : for( iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize )
1787 : {
1788 400 : for( iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize )
1789 : {
1790 : int nTBXSize, nTBYSize;
1791 :
1792 200 : if( !pfnProgress( (nBlocksDone++) / (float) nBlockTotal,
1793 : NULL, pProgressData ) )
1794 : {
1795 : CPLError( CE_Failure, CPLE_UserInterrupt,
1796 0 : "User terminated" );
1797 0 : delete poDS;
1798 0 : CPLFree(pData);
1799 :
1800 : GDALDriver *poHKVDriver =
1801 0 : (GDALDriver *) GDALGetDriverByName( "MFF2" );
1802 0 : poHKVDriver->Delete( pszFilename );
1803 0 : return NULL;
1804 : }
1805 :
1806 200 : nTBXSize = MIN(nBlockXSize,nXSize-iXOffset);
1807 200 : nTBYSize = MIN(nBlockYSize,nYSize-iYOffset);
1808 :
1809 : eErr = poSrcBand->RasterIO( GF_Read,
1810 : iXOffset, iYOffset,
1811 : nTBXSize, nTBYSize,
1812 : pData, nTBXSize, nTBYSize,
1813 200 : eType, 0, 0 );
1814 200 : if( eErr != CE_None )
1815 : {
1816 0 : delete poDS;
1817 0 : CPLFree(pData);
1818 0 : return NULL;
1819 : }
1820 :
1821 : eErr = poDstBand->RasterIO( GF_Write,
1822 : iXOffset, iYOffset,
1823 : nTBXSize, nTBYSize,
1824 : pData, nTBXSize, nTBYSize,
1825 200 : eType, 0, 0 );
1826 :
1827 200 : if( eErr != CE_None )
1828 : {
1829 0 : delete poDS;
1830 0 : CPLFree(pData);
1831 0 : return NULL;
1832 : }
1833 : }
1834 : }
1835 :
1836 20 : CPLFree( pData );
1837 : }
1838 :
1839 : /* -------------------------------------------------------------------- */
1840 : /* Copy georeferencing information, if enough is available. */
1841 : /* Only copy geotransform-style info (won't work for slant range). */
1842 : /* -------------------------------------------------------------------- */
1843 :
1844 10 : double *tempGeoTransform=NULL;
1845 :
1846 10 : tempGeoTransform = (double *) CPLMalloc(6*sizeof(double));
1847 :
1848 20 : if (( poSrcDS->GetGeoTransform( tempGeoTransform ) == CE_None)
1849 10 : && (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0
1850 0 : || tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0
1851 0 : || tempGeoTransform[4] != 0.0 || ABS(tempGeoTransform[5]) != 1.0 ))
1852 : {
1853 :
1854 10 : poDS->SetGCPProjection(poSrcDS->GetProjectionRef());
1855 10 : poDS->SetProjection(poSrcDS->GetProjectionRef());
1856 10 : poDS->SetGeoTransform(tempGeoTransform);
1857 :
1858 10 : CPLFree(tempGeoTransform);
1859 :
1860 : /* georef file will be saved automatically when dataset is deleted */
1861 : /* because SetProjection sets a flag to indicate it's necessary. */
1862 :
1863 : }
1864 : else
1865 : {
1866 0 : CPLFree(tempGeoTransform);
1867 : }
1868 :
1869 : /* Make sure image data gets flushed */
1870 30 : for( iBand = 0; iBand < poDS->GetRasterCount(); iBand++ )
1871 : {
1872 20 : RawRasterBand *poDstBand = (RawRasterBand *) poDS->GetRasterBand( iBand+1 );
1873 20 : poDstBand->FlushCache();
1874 : }
1875 :
1876 :
1877 10 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1878 : {
1879 : CPLError( CE_Failure, CPLE_UserInterrupt,
1880 0 : "User terminated" );
1881 0 : delete poDS;
1882 :
1883 : GDALDriver *poHKVDriver =
1884 0 : (GDALDriver *) GDALGetDriverByName( "MFF2" );
1885 0 : poHKVDriver->Delete( pszFilename );
1886 0 : return NULL;
1887 : }
1888 :
1889 10 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1890 :
1891 10 : return poDS;
1892 : }
1893 :
1894 :
1895 : /************************************************************************/
1896 : /* GDALRegister_HKV() */
1897 : /************************************************************************/
1898 :
1899 582 : void GDALRegister_HKV()
1900 :
1901 : {
1902 : GDALDriver *poDriver;
1903 :
1904 582 : if( GDALGetDriverByName( "MFF2" ) == NULL )
1905 : {
1906 561 : poDriver = new GDALDriver();
1907 :
1908 561 : poDriver->SetDescription( "MFF2" );
1909 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1910 561 : "Vexcel MFF2 (HKV) Raster" );
1911 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1912 561 : "frmt_mff2.html" );
1913 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1914 561 : "Byte Int16 UInt16 Int32 UInt32 CInt16 CInt32 Float32 Float64 CFloat32 CFloat64" );
1915 :
1916 561 : poDriver->pfnOpen = HKVDataset::Open;
1917 561 : poDriver->pfnCreate = HKVDataset::Create;
1918 561 : poDriver->pfnDelete = HKVDataset::Delete;
1919 561 : poDriver->pfnCreateCopy = HKVDataset::CreateCopy;
1920 :
1921 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1922 : }
1923 582 : }
|
