1 : /******************************************************************************
2 : * $Id: hkvdataset.cpp 18183 2009-12-05 01:20:02Z warmerdam $
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 18183 2009-12-05 01:20:02Z warmerdam $");
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, FILE * 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 11 : HKVSpheroidList :: HKVSpheroidList()
79 : {
80 11 : num_spheroids = 58;
81 11 : epsilonR = 0.1;
82 11 : epsilonI = 0.000001;
83 :
84 11 : spheroids[0].SetValuesByEqRadiusAndInvFlattening("airy-1830",6377563.396,299.3249646);
85 11 : spheroids[1].SetValuesByEqRadiusAndInvFlattening("modified-airy",6377340.189,299.3249646);
86 11 : spheroids[2].SetValuesByEqRadiusAndInvFlattening("australian-national",6378160,298.25);
87 11 : spheroids[3].SetValuesByEqRadiusAndInvFlattening("bessel-1841-namibia",6377483.865,299.1528128);
88 11 : spheroids[4].SetValuesByEqRadiusAndInvFlattening("bessel-1841",6377397.155,299.1528128);
89 11 : spheroids[5].SetValuesByEqRadiusAndInvFlattening("clarke-1858",6378294.0,294.297);
90 11 : spheroids[6].SetValuesByEqRadiusAndInvFlattening("clarke-1866",6378206.4,294.9786982);
91 11 : spheroids[7].SetValuesByEqRadiusAndInvFlattening("clarke-1880",6378249.145,293.465);
92 11 : spheroids[8].SetValuesByEqRadiusAndInvFlattening("everest-india-1830",6377276.345,300.8017);
93 11 : spheroids[9].SetValuesByEqRadiusAndInvFlattening("everest-sabah-sarawak",6377298.556,300.8017);
94 11 : spheroids[10].SetValuesByEqRadiusAndInvFlattening("everest-india-1956",6377301.243,300.8017);
95 11 : spheroids[11].SetValuesByEqRadiusAndInvFlattening("everest-malaysia-1969",6377295.664,300.8017);
96 11 : spheroids[12].SetValuesByEqRadiusAndInvFlattening("everest-malay-sing",6377304.063,300.8017);
97 11 : spheroids[13].SetValuesByEqRadiusAndInvFlattening("everest-pakistan",6377309.613,300.8017);
98 11 : spheroids[14].SetValuesByEqRadiusAndInvFlattening("modified-fisher-1960",6378155,298.3);
99 11 : spheroids[15].SetValuesByEqRadiusAndInvFlattening("helmert-1906",6378200,298.3);
100 11 : spheroids[16].SetValuesByEqRadiusAndInvFlattening("hough-1960",6378270,297);
101 11 : spheroids[17].SetValuesByEqRadiusAndInvFlattening("hughes",6378273.0,298.279);
102 11 : spheroids[18].SetValuesByEqRadiusAndInvFlattening("indonesian-1974",6378160,298.247);
103 11 : spheroids[19].SetValuesByEqRadiusAndInvFlattening("international-1924",6378388,297);
104 11 : spheroids[20].SetValuesByEqRadiusAndInvFlattening("iugc-67",6378160.0,298.254);
105 11 : spheroids[21].SetValuesByEqRadiusAndInvFlattening("iugc-75",6378140.0,298.25298);
106 11 : spheroids[22].SetValuesByEqRadiusAndInvFlattening("krassovsky-1940",6378245,298.3);
107 11 : spheroids[23].SetValuesByEqRadiusAndInvFlattening("kaula",6378165.0,292.308);
108 11 : spheroids[24].SetValuesByEqRadiusAndInvFlattening("grs-80",6378137,298.257222101);
109 11 : spheroids[25].SetValuesByEqRadiusAndInvFlattening("south-american-1969",6378160,298.25);
110 11 : spheroids[26].SetValuesByEqRadiusAndInvFlattening("wgs-72",6378135,298.26);
111 11 : spheroids[27].SetValuesByEqRadiusAndInvFlattening("wgs-84",6378137,298.257223563);
112 11 : spheroids[28].SetValuesByEqRadiusAndInvFlattening("ev-wgs-84",6378137.0,298.252841);
113 11 : spheroids[29].SetValuesByEqRadiusAndInvFlattening("ev-bessel",6377397.0,299.1976073);
114 :
115 11 : spheroids[30].SetValuesByEqRadiusAndInvFlattening("airy_1830",6377563.396,299.3249646);
116 11 : spheroids[31].SetValuesByEqRadiusAndInvFlattening("modified_airy",6377340.189,299.3249646);
117 11 : spheroids[32].SetValuesByEqRadiusAndInvFlattening("australian_national",6378160,298.25);
118 11 : spheroids[33].SetValuesByEqRadiusAndInvFlattening("bessel_1841_namibia",6377483.865,299.1528128);
119 11 : spheroids[34].SetValuesByEqRadiusAndInvFlattening("bessel_1841",6377397.155,299.1528128);
120 11 : spheroids[35].SetValuesByEqRadiusAndInvFlattening("clarke_1858",6378294.0,294.297);
121 11 : spheroids[36].SetValuesByEqRadiusAndInvFlattening("clarke_1866",6378206.4,294.9786982);
122 11 : spheroids[37].SetValuesByEqRadiusAndInvFlattening("clarke_1880",6378249.145,293.465);
123 11 : spheroids[38].SetValuesByEqRadiusAndInvFlattening("everest_india_1830",6377276.345,300.8017);
124 11 : spheroids[39].SetValuesByEqRadiusAndInvFlattening("everest_sabah_sarawak",6377298.556,300.8017);
125 11 : spheroids[40].SetValuesByEqRadiusAndInvFlattening("everest_india_1956",6377301.243,300.8017);
126 11 : spheroids[41].SetValuesByEqRadiusAndInvFlattening("everest_malaysia_1969",6377295.664,300.8017);
127 11 : spheroids[42].SetValuesByEqRadiusAndInvFlattening("everest_malay_sing",6377304.063,300.8017);
128 11 : spheroids[43].SetValuesByEqRadiusAndInvFlattening("everest_pakistan",6377309.613,300.8017);
129 11 : spheroids[44].SetValuesByEqRadiusAndInvFlattening("modified_fisher_1960",6378155,298.3);
130 11 : spheroids[45].SetValuesByEqRadiusAndInvFlattening("helmert_1906",6378200,298.3);
131 11 : spheroids[46].SetValuesByEqRadiusAndInvFlattening("hough_1960",6378270,297);
132 11 : spheroids[47].SetValuesByEqRadiusAndInvFlattening("indonesian_1974",6378160,298.247);
133 11 : spheroids[48].SetValuesByEqRadiusAndInvFlattening("international_1924",6378388,297);
134 11 : spheroids[49].SetValuesByEqRadiusAndInvFlattening("iugc_67",6378160.0,298.254);
135 11 : spheroids[50].SetValuesByEqRadiusAndInvFlattening("iugc_75",6378140.0,298.25298);
136 11 : spheroids[51].SetValuesByEqRadiusAndInvFlattening("krassovsky_1940",6378245,298.3);
137 11 : spheroids[52].SetValuesByEqRadiusAndInvFlattening("grs_80",6378137,298.257222101);
138 11 : spheroids[53].SetValuesByEqRadiusAndInvFlattening("south_american_1969",6378160,298.25);
139 11 : spheroids[54].SetValuesByEqRadiusAndInvFlattening("wgs_72",6378135,298.26);
140 11 : spheroids[55].SetValuesByEqRadiusAndInvFlattening("wgs_84",6378137,298.257223563);
141 11 : spheroids[56].SetValuesByEqRadiusAndInvFlattening("ev_wgs_84",6378137.0,298.252841);
142 11 : spheroids[57].SetValuesByEqRadiusAndInvFlattening("ev_bessel",6377397.0,299.1976073);
143 :
144 11 : }
145 :
146 11 : HKVSpheroidList::~HKVSpheroidList()
147 :
148 : {
149 11 : }
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 : FILE *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 41 : HKVRasterBand::HKVRasterBand( HKVDataset *poDS, int nBand, FILE * 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 41 : nLineOffset, eDataType, bNativeOrder, TRUE )
247 :
248 : {
249 41 : this->poDS = poDS;
250 41 : this->nBand = nBand;
251 :
252 41 : nBlockXSize = poDS->GetRasterXSize();
253 41 : nBlockYSize = 1;
254 41 : }
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 82 : HKVRasterBand::~HKVRasterBand()
275 :
276 : {
277 82 : }
278 :
279 : /************************************************************************/
280 : /* ==================================================================== */
281 : /* HKVDataset */
282 : /* ==================================================================== */
283 : /************************************************************************/
284 :
285 : /************************************************************************/
286 : /* HKVDataset() */
287 : /************************************************************************/
288 :
289 21 : HKVDataset::HKVDataset()
290 : {
291 21 : pszPath = NULL;
292 21 : papszAttrib = NULL;
293 21 : papszGeoref = NULL;
294 21 : bGeorefChanged = FALSE;
295 :
296 21 : nGCPCount = 0;
297 21 : pasGCPList = NULL;
298 21 : pszProjection = CPLStrdup("");
299 21 : pszGCPProjection = CPLStrdup("");
300 21 : adfGeoTransform[0] = 0.0;
301 21 : adfGeoTransform[1] = 1.0;
302 21 : adfGeoTransform[2] = 0.0;
303 21 : adfGeoTransform[3] = 0.0;
304 21 : adfGeoTransform[4] = 0.0;
305 21 : adfGeoTransform[5] = 1.0;
306 :
307 21 : bNoDataSet = FALSE;
308 21 : bNoDataChanged = FALSE;
309 :
310 : /* Initialize datasets to new version; change if necessary */
311 21 : MFF2version = (float) 1.1;
312 21 : }
313 :
314 : /************************************************************************/
315 : /* ~HKVDataset() */
316 : /************************************************************************/
317 :
318 42 : HKVDataset::~HKVDataset()
319 :
320 : {
321 21 : FlushCache();
322 21 : if( bGeorefChanged )
323 : {
324 : const char *pszFilename;
325 :
326 10 : pszFilename = CPLFormFilename(pszPath, "georef", NULL );
327 :
328 10 : CSLSave( papszGeoref, pszFilename );
329 : }
330 :
331 21 : 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 21 : if( fpBlob != NULL )
344 21 : VSIFCloseL( fpBlob );
345 :
346 21 : if( nGCPCount > 0 )
347 : {
348 11 : GDALDeinitGCPs( nGCPCount, pasGCPList );
349 11 : CPLFree( pasGCPList );
350 : }
351 :
352 21 : CPLFree( pszProjection );
353 21 : CPLFree( pszGCPProjection );
354 21 : CPLFree( pszPath );
355 21 : CSLDestroy( papszGeoref );
356 21 : CSLDestroy( papszAttrib );
357 42 : }
358 :
359 : /************************************************************************/
360 : /* SetVersion() */
361 : /************************************************************************/
362 :
363 21 : void HKVDataset::SetVersion(float version_number)
364 :
365 : {
366 : //update stored info
367 21 : MFF2version = version_number;
368 21 : }
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 20 : 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 20 : pszFilename = CPLFormFilename( pszFilenameIn, "attrib", NULL );
408 :
409 20 : fp = VSIFOpen( pszFilename, "wt" );
410 20 : 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 20 : fprintf( fp, "channel.enumeration = %d\n", nBands );
418 20 : fprintf( fp, "channel.interleave = { *pixel tile sequential }\n" );
419 20 : fprintf( fp, "extent.cols = %d\n", nXSize );
420 20 : fprintf( fp, "extent.rows = %d\n", nYSize );
421 :
422 20 : 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 2 : "{ *unsigned twos-complement ieee-754 }\n" );
432 2 : break;
433 :
434 : case GDT_CInt16:
435 : case GDT_Int16:
436 : fprintf( fp, "pixel.encoding = "
437 4 : "{ unsigned *twos-complement ieee-754 }\n" );
438 4 : break;
439 :
440 : case GDT_CFloat32:
441 : case GDT_Float32:
442 : fprintf( fp, "pixel.encoding = "
443 4 : "{ unsigned twos-complement *ieee-754 }\n" );
444 : break;
445 :
446 : default:
447 : CPLAssert( FALSE );
448 : }
449 :
450 20 : fprintf( fp, "pixel.size = %d\n", GDALGetDataTypeSize(eType) );
451 20 : if( GDALDataTypeIsComplex( eType ) )
452 4 : fprintf( fp, "pixel.field = { real *complex }\n" );
453 : else
454 16 : fprintf( fp, "pixel.field = { *real complex }\n" );
455 :
456 : #ifdef CPL_MSB
457 : fprintf( fp, "pixel.order = { lsbf *msbf }\n" );
458 : #else
459 20 : fprintf( fp, "pixel.order = { *lsbf msbf }\n" );
460 : #endif
461 :
462 20 : if ( bNoDataSet )
463 0 : fprintf( fp, "pixel.no_data = %f\n", dfNoDataValue );
464 :
465 : /* version information- only create the new style */
466 20 : fprintf( fp, "version = 1.1");
467 :
468 :
469 20 : VSIFClose( fp );
470 20 : 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 10 : CPLErr HKVDataset::GetGeoTransform( double * padfTransform )
489 :
490 : {
491 10 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
492 10 : return( CE_None );
493 : }
494 :
495 : /************************************************************************/
496 : /* SetGeoTransform() */
497 : /************************************************************************/
498 :
499 10 : 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 10 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
512 :
513 : /* Clear previous gcps */
514 10 : if( nGCPCount > 0 )
515 : {
516 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
517 0 : CPLFree( pasGCPList );
518 : }
519 10 : nGCPCount = 0;
520 10 : pasGCPList = NULL;
521 :
522 : /* Return if the identity transform is set */
523 10 : 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 10 : OGRSpatialReference oUTM;
533 10 : OGRSpatialReference oLL;
534 10 : OGRCoordinateTransformation *poTransform = NULL;
535 10 : int bSuccess=TRUE;
536 : char *pszPtemp;
537 : char *pszGCPtemp;
538 :
539 : /* clear old gcps, initialize new list */
540 :
541 10 : if( nGCPCount > 0 )
542 : {
543 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
544 0 : CPLFree( pasGCPList );
545 0 : pasGCPList = NULL;
546 : }
547 10 : nGCPCount = 0;
548 10 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5);
549 :
550 : /* Projection parameter checking will have been done */
551 : /* in SetProjection. */
552 10 : if(( CSLFetchNameValue( papszGeoref, "projection.name" ) != NULL ) &&
553 : ( EQUAL(CSLFetchNameValue( papszGeoref, "projection.name" ),"UTM" )))
554 :
555 : {
556 : /* pass copies of projection info, not originals (pointers */
557 : /* get updated by importFromWkt) */
558 0 : pszPtemp = CPLStrdup(pszProjection);
559 0 : oUTM.importFromWkt(&pszPtemp);
560 0 : (oUTM.GetAttrNode("GEOGCS"))->exportToWkt(&pszGCPtemp);
561 0 : oLL.importFromWkt(&pszGCPtemp);
562 0 : poTransform = OGRCreateCoordinateTransformation( &oUTM, &oLL );
563 0 : if( poTransform == NULL )
564 : {
565 0 : bSuccess = FALSE;
566 0 : CPLErrorReset();
567 : }
568 : }
569 10 : else if ((( CSLFetchNameValue( papszGeoref, "projection.name" ) != NULL ) &&
570 : ( !EQUAL(CSLFetchNameValue( papszGeoref, "projection.name" ),"LL" ))) ||
571 : (CSLFetchNameValue( papszGeoref, "projection.name" ) == NULL ))
572 : {
573 0 : return CE_Failure;
574 : }
575 : /* -------------------------------------------------------------------- */
576 : /* top left */
577 : /* -------------------------------------------------------------------- */
578 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
579 10 : CPLFree( pasGCPList[nGCPCount].pszId );
580 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "top_left" );
581 :
582 10 : if (MFF2version > 1.0)
583 : {
584 10 : temp_lat = padfTransform[3];
585 10 : temp_long = padfTransform[0];
586 10 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
587 10 : pasGCPList[nGCPCount].dfGCPLine = 0.0;
588 : }
589 : else
590 : {
591 0 : temp_lat = padfTransform[3] + 0.5 * padfTransform[4] + 0.5 * padfTransform[5];
592 0 : temp_long = padfTransform[0] + 0.5 * padfTransform[1]+ 0.5 * padfTransform[2];
593 0 : pasGCPList[nGCPCount].dfGCPPixel = 0.5;
594 0 : pasGCPList[nGCPCount].dfGCPLine = 0.5;
595 : }
596 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
597 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
598 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
599 10 : nGCPCount++;
600 :
601 10 : if (poTransform != NULL)
602 : {
603 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
604 0 : bSuccess = FALSE;
605 : }
606 :
607 10 : if (bSuccess)
608 : {
609 10 : sprintf( szValue, "%.10f", temp_lat );
610 : papszGeoref = CSLSetNameValue( papszGeoref, "top_left.latitude",
611 10 : szValue );
612 :
613 10 : sprintf( szValue, "%.10f", temp_long );
614 : papszGeoref = CSLSetNameValue( papszGeoref, "top_left.longitude",
615 10 : szValue );
616 : }
617 :
618 : /* -------------------------------------------------------------------- */
619 : /* top_right */
620 : /* -------------------------------------------------------------------- */
621 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
622 10 : CPLFree( pasGCPList[nGCPCount].pszId );
623 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "top_right" );
624 :
625 10 : if (MFF2version > 1.0)
626 : {
627 10 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4];
628 10 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1];
629 10 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
630 10 : pasGCPList[nGCPCount].dfGCPLine = 0.0;
631 : }
632 : else
633 : {
634 0 : temp_lat = padfTransform[3] + (GetRasterXSize()-0.5) * padfTransform[4] + 0.5 * padfTransform[5];
635 0 : temp_long = padfTransform[0] + (GetRasterXSize()-0.5) * padfTransform[1] + 0.5 * padfTransform[2];
636 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()-0.5;
637 0 : pasGCPList[nGCPCount].dfGCPLine = 0.5;
638 : }
639 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
640 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
641 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
642 10 : nGCPCount++;
643 :
644 10 : if (poTransform != NULL)
645 : {
646 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
647 0 : bSuccess = FALSE;
648 : }
649 :
650 10 : if (bSuccess)
651 : {
652 10 : sprintf( szValue, "%.10f", temp_lat );
653 : papszGeoref = CSLSetNameValue( papszGeoref, "top_right.latitude",
654 10 : szValue );
655 :
656 10 : sprintf( szValue, "%.10f", temp_long );
657 : papszGeoref = CSLSetNameValue( papszGeoref, "top_right.longitude",
658 10 : szValue );
659 : }
660 :
661 : /* -------------------------------------------------------------------- */
662 : /* bottom_left */
663 : /* -------------------------------------------------------------------- */
664 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
665 10 : CPLFree( pasGCPList[nGCPCount].pszId );
666 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "bottom_left" );
667 :
668 10 : if (MFF2version > 1.0)
669 : {
670 10 : temp_lat = padfTransform[3] + GetRasterYSize() * padfTransform[5];
671 10 : temp_long = padfTransform[0] + GetRasterYSize() * padfTransform[2];
672 10 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
673 10 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize();
674 : }
675 : else
676 : {
677 0 : temp_lat = padfTransform[3] + 0.5 * padfTransform[4] + (GetRasterYSize()-0.5) * padfTransform[5];
678 0 : temp_long = padfTransform[0] + 0.5 * padfTransform[1] + (GetRasterYSize()-0.5) * padfTransform[2];
679 0 : pasGCPList[nGCPCount].dfGCPPixel = 0.5;
680 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()-0.5;
681 : }
682 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
683 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
684 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
685 10 : nGCPCount++;
686 :
687 10 : if (poTransform != NULL)
688 : {
689 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
690 0 : bSuccess = FALSE;
691 : }
692 :
693 10 : if (bSuccess)
694 : {
695 10 : sprintf( szValue, "%.10f", temp_lat );
696 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_left.latitude",
697 10 : szValue );
698 :
699 10 : sprintf( szValue, "%.10f", temp_long );
700 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_left.longitude",
701 10 : szValue );
702 : }
703 :
704 : /* -------------------------------------------------------------------- */
705 : /* bottom_right */
706 : /* -------------------------------------------------------------------- */
707 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
708 10 : CPLFree( pasGCPList[nGCPCount].pszId );
709 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "bottom_right" );
710 :
711 10 : if (MFF2version > 1.0)
712 : {
713 20 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] +
714 20 : GetRasterYSize() * padfTransform[5];
715 20 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] +
716 20 : GetRasterYSize() * padfTransform[2];
717 10 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
718 10 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize();
719 :
720 : }
721 : else
722 : {
723 0 : temp_lat = padfTransform[3] + (GetRasterXSize()-0.5) * padfTransform[4] +
724 0 : (GetRasterYSize()-0.5) * padfTransform[5];
725 0 : temp_long = padfTransform[0] + (GetRasterXSize()-0.5) * padfTransform[1] +
726 0 : (GetRasterYSize()-0.5) * padfTransform[2];
727 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()-0.5;
728 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()-0.5;
729 : }
730 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
731 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
732 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
733 10 : nGCPCount++;
734 :
735 10 : if (poTransform != NULL)
736 : {
737 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
738 0 : bSuccess = FALSE;
739 : }
740 :
741 10 : if (bSuccess)
742 : {
743 10 : sprintf( szValue, "%.10f", temp_lat );
744 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_right.latitude",
745 10 : szValue );
746 :
747 10 : sprintf( szValue, "%.10f", temp_long );
748 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_right.longitude",
749 10 : szValue );
750 : }
751 :
752 : /* -------------------------------------------------------------------- */
753 : /* Center */
754 : /* -------------------------------------------------------------------- */
755 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
756 10 : CPLFree( pasGCPList[nGCPCount].pszId );
757 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( "centre" );
758 :
759 10 : if (MFF2version > 1.0)
760 : {
761 20 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] * 0.5 +
762 20 : GetRasterYSize() * padfTransform[5] * 0.5;
763 20 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] * 0.5 +
764 20 : GetRasterYSize() * padfTransform[2] * 0.5;
765 10 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()/2.0;
766 10 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()/2.0;
767 : }
768 : else
769 : {
770 0 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] * 0.5 +
771 0 : GetRasterYSize() * padfTransform[5] * 0.5;
772 0 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] * 0.5 +
773 0 : GetRasterYSize() * padfTransform[2] * 0.5;
774 0 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()/2.0;
775 0 : pasGCPList[nGCPCount].dfGCPLine = GetRasterYSize()/2.0;
776 : }
777 10 : pasGCPList[nGCPCount].dfGCPX = temp_long;
778 10 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
779 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
780 10 : nGCPCount++;
781 :
782 10 : if (poTransform != NULL)
783 : {
784 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
785 0 : bSuccess = FALSE;
786 : }
787 :
788 10 : if (bSuccess)
789 : {
790 10 : sprintf( szValue, "%.10f", temp_lat );
791 : papszGeoref = CSLSetNameValue( papszGeoref, "centre.latitude",
792 10 : szValue );
793 :
794 10 : sprintf( szValue, "%.10f", temp_long );
795 : papszGeoref = CSLSetNameValue( papszGeoref, "centre.longitude",
796 10 : szValue );
797 : }
798 :
799 10 : if (!bSuccess)
800 : {
801 : CPLError(CE_Warning,CPLE_AppDefined,
802 0 : "Warning- error setting header info in SetGeoTransform. Changes may not be saved properly.\n");
803 : }
804 :
805 10 : if (poTransform != NULL)
806 0 : delete poTransform;
807 :
808 :
809 :
810 10 : bGeorefChanged = TRUE;
811 :
812 10 : return( CE_None );
813 : }
814 :
815 10 : CPLErr HKVDataset::SetGCPProjection( const char *pszNewProjection )
816 : {
817 :
818 10 : CPLFree( pszGCPProjection );
819 10 : this->pszGCPProjection = CPLStrdup(pszNewProjection);
820 :
821 10 : return CE_None;
822 : }
823 :
824 : /************************************************************************/
825 : /* SetProjection() */
826 : /* */
827 : /* We provide very limited support for setting the projection. */
828 : /************************************************************************/
829 :
830 10 : CPLErr HKVDataset::SetProjection( const char * pszNewProjection )
831 :
832 : {
833 : HKVSpheroidList *hkvEllipsoids;
834 : double eq_radius, inv_flattening;
835 10 : OGRErr ogrerrorEq=OGRERR_NONE;
836 10 : OGRErr ogrerrorInvf=OGRERR_NONE;
837 10 : OGRErr ogrerrorOl=OGRERR_NONE;
838 :
839 10 : char *spheroid_name = NULL;
840 :
841 : /* This function is used to update a georef file */
842 :
843 :
844 : /* printf( "HKVDataset::SetProjection(%s)\n", pszNewProjection ); */
845 :
846 10 : if( !EQUALN(pszNewProjection,"GEOGCS",6)
847 : && !EQUALN(pszNewProjection,"PROJCS",6)
848 : && !EQUAL(pszNewProjection,"") )
849 : {
850 : CPLError( CE_Failure, CPLE_AppDefined,
851 : "Only OGC WKT Projections supported for writing to HKV.\n"
852 : "%s not supported.",
853 0 : pszNewProjection );
854 :
855 0 : return CE_Failure;
856 : }
857 10 : else if (EQUAL(pszNewProjection,""))
858 : {
859 0 : CPLFree( pszProjection );
860 0 : pszProjection = (char *) CPLStrdup(pszNewProjection);
861 :
862 0 : return CE_None;
863 : }
864 10 : CPLFree( pszProjection );
865 10 : pszProjection = (char *) CPLStrdup(pszNewProjection);
866 :
867 :
868 10 : OGRSpatialReference oSRS(pszNewProjection);
869 :
870 10 : if ((oSRS.GetAttrValue("PROJECTION") != NULL) &&
871 : (EQUAL(oSRS.GetAttrValue("PROJECTION"),SRS_PT_TRANSVERSE_MERCATOR)))
872 : {
873 : char *ol_txt;
874 0 : ol_txt=(char *) CPLMalloc(255);
875 0 : papszGeoref = CSLSetNameValue( papszGeoref, "projection.name", "utm" );
876 0 : sprintf(ol_txt,"%f",oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0,&ogrerrorOl));
877 : papszGeoref = CSLSetNameValue( papszGeoref, "projection.origin_longitude",
878 0 : ol_txt );
879 0 : CPLFree(ol_txt);
880 : }
881 10 : else if ((oSRS.GetAttrValue("PROJECTION") == NULL) && (oSRS.IsGeographic()))
882 : {
883 10 : papszGeoref = CSLSetNameValue( papszGeoref, "projection.name", "LL" );
884 : }
885 : else
886 : {
887 : CPLError( CE_Warning, CPLE_AppDefined,
888 0 : "Unrecognized projection.");
889 0 : return CE_Failure;
890 : }
891 10 : eq_radius = oSRS.GetSemiMajor(&ogrerrorEq);
892 10 : inv_flattening = oSRS.GetInvFlattening(&ogrerrorInvf);
893 20 : if ((ogrerrorEq == OGRERR_NONE) && (ogrerrorInvf == OGRERR_NONE))
894 : {
895 10 : hkvEllipsoids = new HKVSpheroidList;
896 20 : spheroid_name = hkvEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(eq_radius,inv_flattening);
897 10 : if (spheroid_name != NULL)
898 : {
899 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
900 10 : spheroid_name );
901 : }
902 10 : CPLFree(spheroid_name);
903 10 : delete hkvEllipsoids;
904 : }
905 : else
906 : {
907 : /* default to previous behaviour if spheroid not found by */
908 : /* radius and inverse flattening */
909 :
910 0 : if( strstr(pszNewProjection,"Bessel") != NULL )
911 : {
912 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
913 0 : "ev-bessel" );
914 : }
915 : else
916 : {
917 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
918 0 : "ev-wgs-84" );
919 : }
920 : }
921 10 : bGeorefChanged = TRUE;
922 10 : return CE_None;
923 : }
924 :
925 : /************************************************************************/
926 : /* GetGCPCount() */
927 : /************************************************************************/
928 :
929 0 : int HKVDataset::GetGCPCount()
930 :
931 : {
932 0 : return nGCPCount;
933 : }
934 :
935 : /************************************************************************/
936 : /* GetGCPProjection() */
937 : /************************************************************************/
938 :
939 0 : const char *HKVDataset::GetGCPProjection()
940 :
941 : {
942 0 : return pszGCPProjection;
943 : }
944 :
945 : /************************************************************************/
946 : /* GetGCP() */
947 : /************************************************************************/
948 :
949 0 : const GDAL_GCP *HKVDataset::GetGCPs()
950 :
951 : {
952 0 : return pasGCPList;
953 : }
954 :
955 : /************************************************************************/
956 : /* ProcessGeorefGCP() */
957 : /************************************************************************/
958 :
959 5 : void HKVDataset::ProcessGeorefGCP( char **papszGeoref, const char *pszBase,
960 : double dfRasterX, double dfRasterY )
961 :
962 : {
963 : char szFieldName[128];
964 : double dfLat, dfLong;
965 :
966 : /* -------------------------------------------------------------------- */
967 : /* Fetch the GCP from the string list. */
968 : /* -------------------------------------------------------------------- */
969 5 : sprintf( szFieldName, "%s.latitude", pszBase );
970 5 : if( CSLFetchNameValue(papszGeoref, szFieldName) == NULL )
971 0 : return;
972 : else
973 5 : dfLat = atof(CSLFetchNameValue(papszGeoref, szFieldName));
974 :
975 5 : sprintf( szFieldName, "%s.longitude", pszBase );
976 5 : if( CSLFetchNameValue(papszGeoref, szFieldName) == NULL )
977 0 : return;
978 : else
979 5 : dfLong = atof(CSLFetchNameValue(papszGeoref, szFieldName));
980 :
981 : /* -------------------------------------------------------------------- */
982 : /* Add the gcp to the internal list. */
983 : /* -------------------------------------------------------------------- */
984 5 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
985 :
986 5 : CPLFree( pasGCPList[nGCPCount].pszId );
987 :
988 5 : pasGCPList[nGCPCount].pszId = CPLStrdup( pszBase );
989 :
990 5 : pasGCPList[nGCPCount].dfGCPX = dfLong;
991 5 : pasGCPList[nGCPCount].dfGCPY = dfLat;
992 5 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
993 :
994 5 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
995 5 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
996 :
997 5 : nGCPCount++;
998 : }
999 :
1000 : /************************************************************************/
1001 : /* ProcessGeoref() */
1002 : /************************************************************************/
1003 :
1004 1 : void HKVDataset::ProcessGeoref( const char * pszFilename )
1005 :
1006 : {
1007 : int i;
1008 1 : HKVSpheroidList *hkvEllipsoids = NULL;
1009 :
1010 : /* -------------------------------------------------------------------- */
1011 : /* Load the georef file, and boil white space away from around */
1012 : /* the equal sign. */
1013 : /* -------------------------------------------------------------------- */
1014 1 : CSLDestroy( papszGeoref );
1015 1 : papszGeoref = CSLLoad( pszFilename );
1016 1 : if( papszGeoref == NULL )
1017 0 : return;
1018 :
1019 1 : hkvEllipsoids = new HKVSpheroidList;
1020 :
1021 14 : for( i = 0; papszGeoref[i] != NULL; i++ )
1022 : {
1023 13 : int bAfterEqual = FALSE;
1024 : int iSrc, iDst;
1025 13 : char *pszLine = papszGeoref[i];
1026 :
1027 433 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
1028 : {
1029 420 : if( bAfterEqual || pszLine[iSrc] != ' ' )
1030 : {
1031 420 : pszLine[iDst++] = pszLine[iSrc];
1032 : }
1033 :
1034 420 : if( iDst > 0 && pszLine[iDst-1] == '=' )
1035 13 : bAfterEqual = FALSE;
1036 : }
1037 13 : pszLine[iDst] = '\0';
1038 : }
1039 :
1040 : /* -------------------------------------------------------------------- */
1041 : /* Try to get GCPs, in lat/longs . */
1042 : /* -------------------------------------------------------------------- */
1043 1 : nGCPCount = 0;
1044 1 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5);
1045 :
1046 1 : if (MFF2version > 1.0)
1047 : {
1048 : ProcessGeorefGCP( papszGeoref, "top_left",
1049 1 : 0, 0 );
1050 : ProcessGeorefGCP( papszGeoref, "top_right",
1051 1 : GetRasterXSize(), 0 );
1052 : ProcessGeorefGCP( papszGeoref, "bottom_left",
1053 1 : 0, GetRasterYSize() );
1054 : ProcessGeorefGCP( papszGeoref, "bottom_right",
1055 1 : GetRasterXSize(), GetRasterYSize() );
1056 : ProcessGeorefGCP( papszGeoref, "centre",
1057 1 : GetRasterXSize()/2.0, GetRasterYSize()/2.0 );
1058 : }
1059 : else
1060 : {
1061 : ProcessGeorefGCP( papszGeoref, "top_left",
1062 0 : 0.5, 0.5 );
1063 : ProcessGeorefGCP( papszGeoref, "top_right",
1064 0 : GetRasterXSize()-0.5, 0.5 );
1065 : ProcessGeorefGCP( papszGeoref, "bottom_left",
1066 0 : 0.5, GetRasterYSize()-0.5 );
1067 : ProcessGeorefGCP( papszGeoref, "bottom_right",
1068 0 : GetRasterXSize()-0.5, GetRasterYSize()-0.5 );
1069 : ProcessGeorefGCP( papszGeoref, "centre",
1070 0 : GetRasterXSize()/2.0, GetRasterYSize()/2.0 );
1071 : }
1072 :
1073 : /* -------------------------------------------------------------------- */
1074 : /* Do we have a recognised projection? */
1075 : /* -------------------------------------------------------------------- */
1076 : const char *pszProjName, *pszOriginLong, *pszSpheroidName;
1077 : double eq_radius, inv_flattening;
1078 :
1079 : pszProjName = CSLFetchNameValue(papszGeoref,
1080 1 : "projection.name");
1081 : pszOriginLong = CSLFetchNameValue(papszGeoref,
1082 1 : "projection.origin_longitude");
1083 : pszSpheroidName = CSLFetchNameValue(papszGeoref,
1084 1 : "spheroid.name");
1085 :
1086 :
1087 1 : if ((pszSpheroidName != NULL) && (hkvEllipsoids->SpheroidInList(pszSpheroidName)))
1088 : {
1089 1 : eq_radius=hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
1090 1 : inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName);
1091 : }
1092 0 : else if (pszProjName != NULL)
1093 : {
1094 0 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
1095 0 : eq_radius=hkvEllipsoids->GetSpheroidEqRadius("wgs-84");
1096 0 : inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening("wgs-84");
1097 : }
1098 :
1099 2 : if( (pszProjName != NULL) && EQUAL(pszProjName,"utm") && (nGCPCount == 5) )
1100 : {
1101 : /*int nZone = (int)((atof(pszOriginLong)+184.5) / 6.0); */
1102 : int nZone;
1103 :
1104 1 : if (pszOriginLong == NULL)
1105 : {
1106 : /* If origin not specified, assume 0.0 */
1107 : CPLError(CE_Warning,CPLE_AppDefined,
1108 0 : "Warning- no projection origin longitude specified. Assuming 0.0.");
1109 0 : nZone = 31;
1110 : }
1111 : else
1112 1 : nZone = 31 + (int) floor(atof(pszOriginLong)/6.0);
1113 :
1114 1 : OGRSpatialReference oUTM;
1115 1 : OGRSpatialReference oLL;
1116 1 : OGRCoordinateTransformation *poTransform = NULL;
1117 : double dfUtmX[5], dfUtmY[5];
1118 : int gcp_index;
1119 :
1120 1 : int bSuccess = TRUE;
1121 :
1122 1 : if( pasGCPList[4].dfGCPY < 0 )
1123 0 : oUTM.SetUTM( nZone, 0 );
1124 : else
1125 1 : oUTM.SetUTM( nZone, 1 );
1126 :
1127 1 : if (pszOriginLong != NULL)
1128 : {
1129 1 : oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN,atof(pszOriginLong));
1130 1 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));
1131 : }
1132 :
1133 1 : if ((pszSpheroidName == NULL) || (EQUAL(pszSpheroidName,"wgs-84")) ||
1134 : (EQUAL(pszSpheroidName,"wgs_84")))
1135 : {
1136 0 : oUTM.SetWellKnownGeogCS( "WGS84" );
1137 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1138 : }
1139 : else
1140 : {
1141 1 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1142 : {
1143 : oUTM.SetGeogCS( "unknown","unknown",pszSpheroidName,
1144 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1145 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1146 1 : );
1147 : oLL.SetGeogCS( "unknown","unknown",pszSpheroidName,
1148 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1149 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1150 1 : );
1151 : }
1152 : else
1153 : {
1154 0 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
1155 0 : oUTM.SetWellKnownGeogCS( "WGS84" );
1156 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1157 : }
1158 : }
1159 :
1160 1 : poTransform = OGRCreateCoordinateTransformation( &oLL, &oUTM );
1161 1 : if( poTransform == NULL )
1162 : {
1163 0 : CPLErrorReset();
1164 0 : bSuccess = FALSE;
1165 : }
1166 :
1167 6 : for(gcp_index=0;gcp_index<5;gcp_index++)
1168 : {
1169 5 : dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
1170 5 : dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
1171 :
1172 5 : if( bSuccess && !poTransform->Transform( 1, &(dfUtmX[gcp_index]), &(dfUtmY[gcp_index]) ) )
1173 0 : bSuccess = FALSE;
1174 :
1175 : }
1176 :
1177 1 : if( bSuccess )
1178 : {
1179 1 : int transform_ok = FALSE;
1180 :
1181 : /* update GCPS to proper projection */
1182 6 : for(gcp_index=0;gcp_index<5;gcp_index++)
1183 : {
1184 5 : pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
1185 5 : pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
1186 : }
1187 :
1188 1 : CPLFree( pszGCPProjection );
1189 1 : pszGCPProjection = NULL;
1190 1 : oUTM.exportToWkt( &pszGCPProjection );
1191 :
1192 1 : transform_ok = GDALGCPsToGeoTransform(5,pasGCPList,adfGeoTransform,0);
1193 :
1194 1 : CPLFree( pszProjection );
1195 1 : pszProjection = NULL;
1196 1 : if (transform_ok == FALSE)
1197 : {
1198 : /* transform may not be sufficient in all cases (slant range projection) */
1199 0 : adfGeoTransform[0] = 0.0;
1200 0 : adfGeoTransform[1] = 1.0;
1201 0 : adfGeoTransform[2] = 0.0;
1202 0 : adfGeoTransform[3] = 0.0;
1203 0 : adfGeoTransform[4] = 0.0;
1204 0 : adfGeoTransform[5] = 1.0;
1205 0 : pszProjection = CPLStrdup("");
1206 : }
1207 : else
1208 1 : oUTM.exportToWkt( &pszProjection );
1209 :
1210 : }
1211 :
1212 1 : if( poTransform != NULL )
1213 1 : delete poTransform;
1214 : }
1215 0 : else if ((pszProjName != NULL) && (nGCPCount == 5))
1216 : {
1217 0 : OGRSpatialReference oLL;
1218 0 : int transform_ok = FALSE;
1219 :
1220 :
1221 0 : if (pszOriginLong != NULL)
1222 : {
1223 0 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));
1224 : }
1225 :
1226 0 : if ((pszSpheroidName == NULL) || (EQUAL(pszSpheroidName,"wgs-84")) ||
1227 : (EQUAL(pszSpheroidName,"wgs_84")))
1228 : {
1229 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1230 : }
1231 : else
1232 : {
1233 0 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1234 : {
1235 : oLL.SetGeogCS( "","",pszSpheroidName,
1236 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1237 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1238 0 : );
1239 : }
1240 : else
1241 : {
1242 0 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
1243 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1244 : }
1245 : }
1246 :
1247 0 : transform_ok = GDALGCPsToGeoTransform(5,pasGCPList,adfGeoTransform,0);
1248 :
1249 0 : CPLFree( pszProjection );
1250 0 : pszProjection = NULL;
1251 :
1252 0 : if (transform_ok == FALSE)
1253 : {
1254 0 : adfGeoTransform[0] = 0.0;
1255 0 : adfGeoTransform[1] = 1.0;
1256 0 : adfGeoTransform[2] = 0.0;
1257 0 : adfGeoTransform[3] = 0.0;
1258 0 : adfGeoTransform[4] = 0.0;
1259 0 : adfGeoTransform[5] = 1.0;
1260 : }
1261 : else
1262 : {
1263 0 : oLL.exportToWkt( &pszProjection );
1264 : }
1265 :
1266 0 : CPLFree( pszGCPProjection );
1267 0 : pszGCPProjection = NULL;
1268 0 : oLL.exportToWkt( &pszGCPProjection );
1269 :
1270 : }
1271 :
1272 1 : delete hkvEllipsoids;
1273 : }
1274 :
1275 : /************************************************************************/
1276 : /* Open() */
1277 : /************************************************************************/
1278 :
1279 8580 : GDALDataset *HKVDataset::Open( GDALOpenInfo * poOpenInfo )
1280 :
1281 : {
1282 8580 : int i, bNoDataSet = FALSE;
1283 8580 : double dfNoDataValue = 0.0;
1284 : char **papszAttrib;
1285 : const char *pszFilename, *pszValue;
1286 : VSIStatBuf sStat;
1287 :
1288 : /* -------------------------------------------------------------------- */
1289 : /* We assume the dataset is passed as a directory. Check for */
1290 : /* an attrib and blob file as a minimum. */
1291 : /* -------------------------------------------------------------------- */
1292 8580 : if( !poOpenInfo->bIsDirectory )
1293 8558 : return NULL;
1294 :
1295 22 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "image_data", NULL);
1296 22 : if( VSIStat(pszFilename,&sStat) != 0 )
1297 1 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "blob", NULL );
1298 22 : if( VSIStat(pszFilename,&sStat) != 0 )
1299 1 : return NULL;
1300 :
1301 21 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "attrib", NULL );
1302 21 : if( VSIStat(pszFilename,&sStat) != 0 )
1303 0 : return NULL;
1304 :
1305 : /* -------------------------------------------------------------------- */
1306 : /* Load the attrib file, and boil white space away from around */
1307 : /* the equal sign. */
1308 : /* -------------------------------------------------------------------- */
1309 21 : papszAttrib = CSLLoad( pszFilename );
1310 21 : if( papszAttrib == NULL )
1311 0 : return NULL;
1312 :
1313 210 : for( i = 0; papszAttrib[i] != NULL; i++ )
1314 : {
1315 189 : int bAfterEqual = FALSE;
1316 : int iSrc, iDst;
1317 189 : char *pszLine = papszAttrib[i];
1318 :
1319 5322 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
1320 : {
1321 5133 : if( bAfterEqual || pszLine[iSrc] != ' ' )
1322 : {
1323 4461 : pszLine[iDst++] = pszLine[iSrc];
1324 : }
1325 :
1326 5133 : if( iDst > 0 && pszLine[iDst-1] == '=' )
1327 378 : bAfterEqual = FALSE;
1328 : }
1329 189 : pszLine[iDst] = '\0';
1330 : }
1331 :
1332 : /* -------------------------------------------------------------------- */
1333 : /* Create a corresponding GDALDataset. */
1334 : /* -------------------------------------------------------------------- */
1335 : HKVDataset *poDS;
1336 :
1337 21 : poDS = new HKVDataset();
1338 :
1339 21 : poDS->pszPath = CPLStrdup( poOpenInfo->pszFilename );
1340 21 : poDS->papszAttrib = papszAttrib;
1341 :
1342 21 : poDS->eAccess = poOpenInfo->eAccess;
1343 :
1344 : /* -------------------------------------------------------------------- */
1345 : /* Set some dataset wide information. */
1346 : /* -------------------------------------------------------------------- */
1347 : int bNative, bComplex;
1348 21 : int nRawBands = 0;
1349 :
1350 42 : if( CSLFetchNameValue( papszAttrib, "extent.cols" ) == NULL
1351 : || CSLFetchNameValue( papszAttrib, "extent.rows" ) == NULL )
1352 0 : return NULL;
1353 :
1354 21 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib,"extent.cols"));
1355 21 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib,"extent.rows"));
1356 :
1357 21 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
1358 : {
1359 0 : delete poDS;
1360 0 : return NULL;
1361 : }
1362 :
1363 21 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.order");
1364 21 : if( pszValue == NULL )
1365 0 : bNative = TRUE;
1366 : else
1367 : {
1368 : #ifdef CPL_MSB
1369 : bNative = (strstr(pszValue,"*msbf") != NULL);
1370 : #else
1371 21 : bNative = (strstr(pszValue,"*lsbf") != NULL);
1372 : #endif
1373 : }
1374 :
1375 21 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.no_data");
1376 21 : if( pszValue != NULL )
1377 : {
1378 0 : bNoDataSet = TRUE;
1379 0 : dfNoDataValue = atof(pszValue);
1380 : }
1381 :
1382 21 : pszValue = CSLFetchNameValue(papszAttrib,"channel.enumeration");
1383 21 : if( pszValue != NULL )
1384 21 : nRawBands = atoi(pszValue);
1385 : else
1386 0 : nRawBands = 1;
1387 :
1388 21 : if (!GDALCheckBandCount(nRawBands, TRUE))
1389 : {
1390 0 : delete poDS;
1391 0 : return NULL;
1392 : }
1393 :
1394 21 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.field");
1395 25 : if( pszValue != NULL && strstr(pszValue,"*complex") != NULL )
1396 4 : bComplex = TRUE;
1397 : else
1398 17 : bComplex = FALSE;
1399 :
1400 : /* Get the version number, if present (if not, assume old version. */
1401 : /* Versions differ in their interpretation of corner coordinates. */
1402 :
1403 21 : if (CSLFetchNameValue( papszAttrib, "version" ) != NULL)
1404 : poDS->SetVersion((float)
1405 21 : atof(CSLFetchNameValue(papszAttrib, "version")));
1406 : else
1407 0 : poDS->SetVersion(1.0);
1408 :
1409 : /* -------------------------------------------------------------------- */
1410 : /* Figure out the datatype */
1411 : /* -------------------------------------------------------------------- */
1412 : const char * pszEncoding;
1413 21 : int nSize = 1;
1414 : int nPseudoBands;
1415 : GDALDataType eType;
1416 :
1417 21 : pszEncoding = CSLFetchNameValue(papszAttrib,"pixel.encoding");
1418 21 : if( pszEncoding == NULL )
1419 0 : pszEncoding = "{ *unsigned }";
1420 :
1421 21 : if( CSLFetchNameValue(papszAttrib,"pixel.size") != NULL )
1422 21 : nSize = atoi(CSLFetchNameValue(papszAttrib,"pixel.size"))/8;
1423 :
1424 21 : if( bComplex )
1425 4 : nPseudoBands = 2;
1426 : else
1427 17 : nPseudoBands = 1;
1428 :
1429 21 : if( nSize == 1 )
1430 11 : eType = GDT_Byte;
1431 12 : else if( nSize == 2 && strstr(pszEncoding,"*unsigned") != NULL )
1432 2 : eType = GDT_UInt16;
1433 10 : else if( nSize == 4 && bComplex )
1434 2 : eType = GDT_CInt16;
1435 6 : else if( nSize == 2 )
1436 2 : eType = GDT_Int16;
1437 4 : else if( nSize == 4 && strstr(pszEncoding,"*unsigned") != NULL )
1438 0 : eType = GDT_UInt32;
1439 4 : else if( nSize == 8 && strstr(pszEncoding,"*two") != NULL && bComplex )
1440 0 : eType = GDT_CInt32;
1441 4 : else if( nSize == 4 && strstr(pszEncoding,"*two") != NULL )
1442 0 : eType = GDT_Int32;
1443 6 : else if( nSize == 8 && bComplex )
1444 2 : eType = GDT_CFloat32;
1445 2 : else if( nSize == 4 )
1446 2 : eType = GDT_Float32;
1447 0 : else if( nSize == 16 && bComplex )
1448 0 : eType = GDT_CFloat64;
1449 0 : else if( nSize == 8 )
1450 0 : eType = GDT_Float64;
1451 : else
1452 : {
1453 : CPLError( CE_Failure, CPLE_AppDefined,
1454 : "Unsupported pixel data type in %s.\n"
1455 : "pixel.size=%d pixel.encoding=%s\n",
1456 0 : poDS->pszPath, nSize, pszEncoding );
1457 0 : delete poDS;
1458 0 : return NULL;
1459 : }
1460 :
1461 : /* -------------------------------------------------------------------- */
1462 : /* Open the blob file. */
1463 : /* -------------------------------------------------------------------- */
1464 21 : pszFilename = CPLFormFilename(poDS->pszPath, "image_data", NULL );
1465 21 : if( VSIStat(pszFilename,&sStat) != 0 )
1466 0 : pszFilename = CPLFormFilename(poDS->pszPath, "blob", NULL );
1467 21 : if( poOpenInfo->eAccess == GA_ReadOnly )
1468 : {
1469 1 : poDS->fpBlob = VSIFOpenL( pszFilename, "rb" );
1470 1 : if( poDS->fpBlob == NULL )
1471 : {
1472 : CPLError( CE_Failure, CPLE_OpenFailed,
1473 : "Unable to open file %s for read access.\n",
1474 0 : pszFilename );
1475 0 : delete poDS;
1476 0 : return NULL;
1477 : }
1478 : }
1479 : else
1480 : {
1481 20 : poDS->fpBlob = VSIFOpenL( pszFilename, "rb+" );
1482 20 : if( poDS->fpBlob == NULL )
1483 : {
1484 : CPLError( CE_Failure, CPLE_OpenFailed,
1485 : "Unable to open file %s for update access.\n",
1486 0 : pszFilename );
1487 0 : delete poDS;
1488 0 : return NULL;
1489 : }
1490 : }
1491 :
1492 : /* -------------------------------------------------------------------- */
1493 : /* Build the overview filename, as blob file = "_ovr". */
1494 : /* -------------------------------------------------------------------- */
1495 21 : char *pszOvrFilename = (char *) CPLMalloc(strlen(pszFilename)+5);
1496 :
1497 21 : sprintf( pszOvrFilename, "%s_ovr", pszFilename );
1498 :
1499 : /* -------------------------------------------------------------------- */
1500 : /* Define the bands. */
1501 : /* -------------------------------------------------------------------- */
1502 : int nPixelOffset, nLineOffset, nOffset;
1503 :
1504 21 : nPixelOffset = nRawBands * nSize;
1505 21 : nLineOffset = nPixelOffset * poDS->GetRasterXSize();
1506 21 : nOffset = 0;
1507 :
1508 62 : for( int iRawBand=0; iRawBand < nRawBands; iRawBand++ )
1509 : {
1510 : HKVRasterBand *poBand;
1511 :
1512 : poBand =
1513 : new HKVRasterBand( poDS, poDS->GetRasterCount()+1, poDS->fpBlob,
1514 : nOffset, nPixelOffset, nLineOffset,
1515 41 : eType, bNative );
1516 41 : poDS->SetBand( poDS->GetRasterCount()+1, poBand );
1517 41 : nOffset += GDALGetDataTypeSize( eType ) / 8;
1518 :
1519 41 : if( bNoDataSet )
1520 0 : poBand->SetNoDataValue( dfNoDataValue );
1521 : }
1522 :
1523 21 : poDS->eRasterType = eType;
1524 :
1525 : /* -------------------------------------------------------------------- */
1526 : /* Process the georef file if there is one. */
1527 : /* -------------------------------------------------------------------- */
1528 21 : pszFilename = CPLFormFilename(poDS->pszPath, "georef", NULL );
1529 21 : if( VSIStat(pszFilename,&sStat) == 0 )
1530 1 : poDS->ProcessGeoref(pszFilename);
1531 :
1532 : /* -------------------------------------------------------------------- */
1533 : /* Initialize any PAM information. */
1534 : /* -------------------------------------------------------------------- */
1535 21 : poDS->SetDescription( pszOvrFilename );
1536 21 : poDS->TryLoadXML();
1537 :
1538 : /* -------------------------------------------------------------------- */
1539 : /* Handle overviews. */
1540 : /* -------------------------------------------------------------------- */
1541 21 : poDS->oOvManager.Initialize( poDS, pszOvrFilename, NULL, TRUE );
1542 :
1543 21 : CPLFree( pszOvrFilename );
1544 :
1545 21 : return( poDS );
1546 : }
1547 :
1548 : /************************************************************************/
1549 : /* Create() */
1550 : /************************************************************************/
1551 :
1552 31 : GDALDataset *HKVDataset::Create( const char * pszFilenameIn,
1553 : int nXSize, int nYSize, int nBands,
1554 : GDALDataType eType,
1555 : char ** /* papszParmList */ )
1556 :
1557 : {
1558 : /* -------------------------------------------------------------------- */
1559 : /* Verify input options. */
1560 : /* -------------------------------------------------------------------- */
1561 31 : if (nBands <= 0)
1562 : {
1563 : CPLError( CE_Failure, CPLE_NotSupported,
1564 1 : "HKV driver does not support %d bands.\n", nBands);
1565 1 : return NULL;
1566 : }
1567 :
1568 30 : if( eType != GDT_Byte
1569 : && eType != GDT_UInt16 && eType != GDT_Int16
1570 : && eType != GDT_CInt16 && eType != GDT_Float32
1571 : && eType != GDT_CFloat32 )
1572 : {
1573 : CPLError( CE_Failure, CPLE_AppDefined,
1574 : "Attempt to create HKV file with currently unsupported\n"
1575 : "data type (%s).\n",
1576 10 : GDALGetDataTypeName(eType) );
1577 :
1578 10 : return NULL;
1579 : }
1580 :
1581 : /* -------------------------------------------------------------------- */
1582 : /* Establish the name of the directory we will be creating the */
1583 : /* new HKV directory in. Verify that this is a directory. */
1584 : /* -------------------------------------------------------------------- */
1585 : char *pszBaseDir;
1586 : VSIStatBuf sStat;
1587 :
1588 20 : if( strlen(CPLGetPath(pszFilenameIn)) == 0 )
1589 0 : pszBaseDir = CPLStrdup(".");
1590 : else
1591 20 : pszBaseDir = CPLStrdup(CPLGetPath(pszFilenameIn));
1592 :
1593 20 : if( CPLStat( pszBaseDir, &sStat ) != 0 || !VSI_ISDIR( sStat.st_mode ) )
1594 : {
1595 : CPLError( CE_Failure, CPLE_AppDefined,
1596 : "Attempt to create HKV dataset under %s,\n"
1597 : "but this is not a valid directory.\n",
1598 0 : pszBaseDir);
1599 0 : CPLFree( pszBaseDir );
1600 0 : return NULL;
1601 : }
1602 :
1603 20 : if( VSIMkdir( pszFilenameIn, 0755 ) != 0 )
1604 : {
1605 : CPLError( CE_Failure, CPLE_AppDefined,
1606 : "Unable to create directory %s.\n",
1607 0 : pszFilenameIn );
1608 0 : return NULL;
1609 : }
1610 :
1611 20 : CPLFree( pszBaseDir );
1612 :
1613 : /* -------------------------------------------------------------------- */
1614 : /* Create the header file. */
1615 : /* -------------------------------------------------------------------- */
1616 : CPLErr CEHeaderCreated;
1617 :
1618 : CEHeaderCreated = SaveHKVAttribFile( pszFilenameIn, nXSize, nYSize,
1619 20 : nBands, eType, FALSE, 0.0 );
1620 :
1621 20 : 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 20 : pszFilename = CPLFormFilename( pszFilenameIn, "image_data", NULL );
1632 20 : fp = VSIFOpen( pszFilename, "wb" );
1633 20 : 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 20 : VSIFWrite( (void*)"", 1, 1, fp );
1641 20 : VSIFClose( fp );
1642 :
1643 : /* -------------------------------------------------------------------- */
1644 : /* Open the dataset normally. */
1645 : /* -------------------------------------------------------------------- */
1646 20 : 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 16 : 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 16 : int nBands = poSrcDS->GetRasterCount();
1724 16 : 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 15 : eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1732 :
1733 15 : 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 25 : for( iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++ )
1739 : {
1740 10 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1741 10 : eType = GDALDataTypeUnion( eType, poBand->GetRasterDataType() );
1742 : }
1743 :
1744 : poDS = (HKVDataset *) Create( pszFilename,
1745 : poSrcDS->GetRasterXSize(),
1746 : poSrcDS->GetRasterYSize(),
1747 : poSrcDS->GetRasterCount(),
1748 15 : eType, papszOptions );
1749 :
1750 : /* Check that Create worked- return Null if it didn't */
1751 15 : if (poDS == NULL)
1752 5 : 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 :
1799 : GDALDriver *poHKVDriver =
1800 0 : (GDALDriver *) GDALGetDriverByName( "MFF2" );
1801 0 : poHKVDriver->Delete( pszFilename );
1802 0 : return NULL;
1803 : }
1804 :
1805 200 : nTBXSize = MIN(nBlockXSize,nXSize-iXOffset);
1806 200 : nTBYSize = MIN(nBlockYSize,nYSize-iYOffset);
1807 :
1808 : eErr = poSrcBand->RasterIO( GF_Read,
1809 : iXOffset, iYOffset,
1810 : nTBXSize, nTBYSize,
1811 : pData, nTBXSize, nTBYSize,
1812 200 : eType, 0, 0 );
1813 200 : if( eErr != CE_None )
1814 : {
1815 0 : return NULL;
1816 : }
1817 :
1818 : eErr = poDstBand->RasterIO( GF_Write,
1819 : iXOffset, iYOffset,
1820 : nTBXSize, nTBYSize,
1821 : pData, nTBXSize, nTBYSize,
1822 200 : eType, 0, 0 );
1823 :
1824 200 : if( eErr != CE_None )
1825 : {
1826 0 : return NULL;
1827 : }
1828 : }
1829 : }
1830 :
1831 20 : CPLFree( pData );
1832 : }
1833 :
1834 : /* -------------------------------------------------------------------- */
1835 : /* Copy georeferencing information, if enough is available. */
1836 : /* Only copy geotransform-style info (won't work for slant range). */
1837 : /* -------------------------------------------------------------------- */
1838 :
1839 10 : double *tempGeoTransform=NULL;
1840 :
1841 10 : tempGeoTransform = (double *) CPLMalloc(6*sizeof(double));
1842 :
1843 20 : if (( poSrcDS->GetGeoTransform( tempGeoTransform ) == CE_None)
1844 10 : && (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0
1845 0 : || tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0
1846 0 : || tempGeoTransform[4] != 0.0 || ABS(tempGeoTransform[5]) != 1.0 ))
1847 : {
1848 :
1849 10 : poDS->SetGCPProjection(poSrcDS->GetProjectionRef());
1850 10 : poDS->SetProjection(poSrcDS->GetProjectionRef());
1851 10 : poDS->SetGeoTransform(tempGeoTransform);
1852 :
1853 10 : CPLFree(tempGeoTransform);
1854 :
1855 : /* georef file will be saved automatically when dataset is deleted */
1856 : /* because SetProjection sets a flag to indicate it's necessary. */
1857 :
1858 : }
1859 : else
1860 : {
1861 0 : CPLFree(tempGeoTransform);
1862 : }
1863 :
1864 : /* Make sure image data gets flushed */
1865 30 : for( iBand = 0; iBand < poDS->GetRasterCount(); iBand++ )
1866 : {
1867 20 : RawRasterBand *poDstBand = (RawRasterBand *) poDS->GetRasterBand( iBand+1 );
1868 20 : poDstBand->FlushCache();
1869 : }
1870 :
1871 :
1872 10 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1873 : {
1874 : CPLError( CE_Failure, CPLE_UserInterrupt,
1875 0 : "User terminated" );
1876 0 : delete poDS;
1877 :
1878 : GDALDriver *poHKVDriver =
1879 0 : (GDALDriver *) GDALGetDriverByName( "MFF2" );
1880 0 : poHKVDriver->Delete( pszFilename );
1881 0 : return NULL;
1882 : }
1883 :
1884 10 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1885 :
1886 10 : return poDS;
1887 : }
1888 :
1889 :
1890 : /************************************************************************/
1891 : /* GDALRegister_HKV() */
1892 : /************************************************************************/
1893 :
1894 338 : void GDALRegister_HKV()
1895 :
1896 : {
1897 : GDALDriver *poDriver;
1898 :
1899 338 : if( GDALGetDriverByName( "MFF2" ) == NULL )
1900 : {
1901 336 : poDriver = new GDALDriver();
1902 :
1903 336 : poDriver->SetDescription( "MFF2" );
1904 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1905 336 : "Vexcel MFF2 (HKV) Raster" );
1906 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1907 336 : "frmt_mff2.html" );
1908 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1909 336 : "Byte Int16 UInt16 Int32 UInt32 CInt16 CInt32 Float32 Float64 CFloat32 CFloat64" );
1910 :
1911 336 : poDriver->pfnOpen = HKVDataset::Open;
1912 336 : poDriver->pfnCreate = HKVDataset::Create;
1913 336 : poDriver->pfnDelete = HKVDataset::Delete;
1914 336 : poDriver->pfnCreateCopy = HKVDataset::CreateCopy;
1915 :
1916 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
1917 : }
1918 338 : }
|
