1 : /******************************************************************************
2 : * $Id: hkvdataset.cpp 20996 2010-10-28 18:38:15Z 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 20996 2010-10-28 18:38:15Z 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 82 : HKVSpheroidList :: HKVSpheroidList()
79 : {
80 82 : num_spheroids = 58;
81 82 : epsilonR = 0.1;
82 82 : epsilonI = 0.000001;
83 :
84 82 : spheroids[0].SetValuesByEqRadiusAndInvFlattening("airy-1830",6377563.396,299.3249646);
85 82 : spheroids[1].SetValuesByEqRadiusAndInvFlattening("modified-airy",6377340.189,299.3249646);
86 82 : spheroids[2].SetValuesByEqRadiusAndInvFlattening("australian-national",6378160,298.25);
87 82 : spheroids[3].SetValuesByEqRadiusAndInvFlattening("bessel-1841-namibia",6377483.865,299.1528128);
88 82 : spheroids[4].SetValuesByEqRadiusAndInvFlattening("bessel-1841",6377397.155,299.1528128);
89 82 : spheroids[5].SetValuesByEqRadiusAndInvFlattening("clarke-1858",6378294.0,294.297);
90 82 : spheroids[6].SetValuesByEqRadiusAndInvFlattening("clarke-1866",6378206.4,294.9786982);
91 82 : spheroids[7].SetValuesByEqRadiusAndInvFlattening("clarke-1880",6378249.145,293.465);
92 82 : spheroids[8].SetValuesByEqRadiusAndInvFlattening("everest-india-1830",6377276.345,300.8017);
93 82 : spheroids[9].SetValuesByEqRadiusAndInvFlattening("everest-sabah-sarawak",6377298.556,300.8017);
94 82 : spheroids[10].SetValuesByEqRadiusAndInvFlattening("everest-india-1956",6377301.243,300.8017);
95 82 : spheroids[11].SetValuesByEqRadiusAndInvFlattening("everest-malaysia-1969",6377295.664,300.8017);
96 82 : spheroids[12].SetValuesByEqRadiusAndInvFlattening("everest-malay-sing",6377304.063,300.8017);
97 82 : spheroids[13].SetValuesByEqRadiusAndInvFlattening("everest-pakistan",6377309.613,300.8017);
98 82 : spheroids[14].SetValuesByEqRadiusAndInvFlattening("modified-fisher-1960",6378155,298.3);
99 82 : spheroids[15].SetValuesByEqRadiusAndInvFlattening("helmert-1906",6378200,298.3);
100 82 : spheroids[16].SetValuesByEqRadiusAndInvFlattening("hough-1960",6378270,297);
101 82 : spheroids[17].SetValuesByEqRadiusAndInvFlattening("hughes",6378273.0,298.279);
102 82 : spheroids[18].SetValuesByEqRadiusAndInvFlattening("indonesian-1974",6378160,298.247);
103 82 : spheroids[19].SetValuesByEqRadiusAndInvFlattening("international-1924",6378388,297);
104 82 : spheroids[20].SetValuesByEqRadiusAndInvFlattening("iugc-67",6378160.0,298.254);
105 82 : spheroids[21].SetValuesByEqRadiusAndInvFlattening("iugc-75",6378140.0,298.25298);
106 82 : spheroids[22].SetValuesByEqRadiusAndInvFlattening("krassovsky-1940",6378245,298.3);
107 82 : spheroids[23].SetValuesByEqRadiusAndInvFlattening("kaula",6378165.0,292.308);
108 82 : spheroids[24].SetValuesByEqRadiusAndInvFlattening("grs-80",6378137,298.257222101);
109 82 : spheroids[25].SetValuesByEqRadiusAndInvFlattening("south-american-1969",6378160,298.25);
110 82 : spheroids[26].SetValuesByEqRadiusAndInvFlattening("wgs-72",6378135,298.26);
111 82 : spheroids[27].SetValuesByEqRadiusAndInvFlattening("wgs-84",6378137,298.257223563);
112 82 : spheroids[28].SetValuesByEqRadiusAndInvFlattening("ev-wgs-84",6378137.0,298.252841);
113 82 : spheroids[29].SetValuesByEqRadiusAndInvFlattening("ev-bessel",6377397.0,299.1976073);
114 :
115 82 : spheroids[30].SetValuesByEqRadiusAndInvFlattening("airy_1830",6377563.396,299.3249646);
116 82 : spheroids[31].SetValuesByEqRadiusAndInvFlattening("modified_airy",6377340.189,299.3249646);
117 82 : spheroids[32].SetValuesByEqRadiusAndInvFlattening("australian_national",6378160,298.25);
118 82 : spheroids[33].SetValuesByEqRadiusAndInvFlattening("bessel_1841_namibia",6377483.865,299.1528128);
119 82 : spheroids[34].SetValuesByEqRadiusAndInvFlattening("bessel_1841",6377397.155,299.1528128);
120 82 : spheroids[35].SetValuesByEqRadiusAndInvFlattening("clarke_1858",6378294.0,294.297);
121 82 : spheroids[36].SetValuesByEqRadiusAndInvFlattening("clarke_1866",6378206.4,294.9786982);
122 82 : spheroids[37].SetValuesByEqRadiusAndInvFlattening("clarke_1880",6378249.145,293.465);
123 82 : spheroids[38].SetValuesByEqRadiusAndInvFlattening("everest_india_1830",6377276.345,300.8017);
124 82 : spheroids[39].SetValuesByEqRadiusAndInvFlattening("everest_sabah_sarawak",6377298.556,300.8017);
125 82 : spheroids[40].SetValuesByEqRadiusAndInvFlattening("everest_india_1956",6377301.243,300.8017);
126 82 : spheroids[41].SetValuesByEqRadiusAndInvFlattening("everest_malaysia_1969",6377295.664,300.8017);
127 82 : spheroids[42].SetValuesByEqRadiusAndInvFlattening("everest_malay_sing",6377304.063,300.8017);
128 82 : spheroids[43].SetValuesByEqRadiusAndInvFlattening("everest_pakistan",6377309.613,300.8017);
129 82 : spheroids[44].SetValuesByEqRadiusAndInvFlattening("modified_fisher_1960",6378155,298.3);
130 82 : spheroids[45].SetValuesByEqRadiusAndInvFlattening("helmert_1906",6378200,298.3);
131 82 : spheroids[46].SetValuesByEqRadiusAndInvFlattening("hough_1960",6378270,297);
132 82 : spheroids[47].SetValuesByEqRadiusAndInvFlattening("indonesian_1974",6378160,298.247);
133 82 : spheroids[48].SetValuesByEqRadiusAndInvFlattening("international_1924",6378388,297);
134 82 : spheroids[49].SetValuesByEqRadiusAndInvFlattening("iugc_67",6378160.0,298.254);
135 82 : spheroids[50].SetValuesByEqRadiusAndInvFlattening("iugc_75",6378140.0,298.25298);
136 82 : spheroids[51].SetValuesByEqRadiusAndInvFlattening("krassovsky_1940",6378245,298.3);
137 82 : spheroids[52].SetValuesByEqRadiusAndInvFlattening("grs_80",6378137,298.257222101);
138 82 : spheroids[53].SetValuesByEqRadiusAndInvFlattening("south_american_1969",6378160,298.25);
139 82 : spheroids[54].SetValuesByEqRadiusAndInvFlattening("wgs_72",6378135,298.26);
140 82 : spheroids[55].SetValuesByEqRadiusAndInvFlattening("wgs_84",6378137,298.257223563);
141 82 : spheroids[56].SetValuesByEqRadiusAndInvFlattening("ev_wgs_84",6378137.0,298.252841);
142 82 : spheroids[57].SetValuesByEqRadiusAndInvFlattening("ev_bessel",6377397.0,299.1976073);
143 :
144 82 : }
145 :
146 82 : HKVSpheroidList::~HKVSpheroidList()
147 :
148 : {
149 82 : }
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 182 : 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 182 : nLineOffset, eDataType, bNativeOrder, TRUE )
247 :
248 : {
249 182 : this->poDS = poDS;
250 182 : this->nBand = nBand;
251 :
252 182 : nBlockXSize = poDS->GetRasterXSize();
253 182 : nBlockYSize = 1;
254 182 : }
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 182 : HKVRasterBand::~HKVRasterBand()
275 :
276 : {
277 182 : }
278 :
279 : /************************************************************************/
280 : /* ==================================================================== */
281 : /* HKVDataset */
282 : /* ==================================================================== */
283 : /************************************************************************/
284 :
285 : /************************************************************************/
286 : /* HKVDataset() */
287 : /************************************************************************/
288 :
289 82 : HKVDataset::HKVDataset()
290 : {
291 82 : pszPath = NULL;
292 82 : papszAttrib = NULL;
293 82 : papszGeoref = NULL;
294 82 : bGeorefChanged = FALSE;
295 :
296 82 : nGCPCount = 0;
297 82 : pasGCPList = NULL;
298 82 : pszProjection = CPLStrdup("");
299 82 : pszGCPProjection = CPLStrdup("");
300 82 : adfGeoTransform[0] = 0.0;
301 82 : adfGeoTransform[1] = 1.0;
302 82 : adfGeoTransform[2] = 0.0;
303 82 : adfGeoTransform[3] = 0.0;
304 82 : adfGeoTransform[4] = 0.0;
305 82 : adfGeoTransform[5] = 1.0;
306 :
307 82 : bNoDataSet = FALSE;
308 82 : bNoDataChanged = FALSE;
309 :
310 : /* Initialize datasets to new version; change if necessary */
311 82 : MFF2version = (float) 1.1;
312 82 : }
313 :
314 : /************************************************************************/
315 : /* ~HKVDataset() */
316 : /************************************************************************/
317 :
318 82 : HKVDataset::~HKVDataset()
319 :
320 : {
321 82 : FlushCache();
322 82 : if( bGeorefChanged )
323 : {
324 : const char *pszFilename;
325 :
326 50 : pszFilename = CPLFormFilename(pszPath, "georef", NULL );
327 :
328 50 : CSLSave( papszGeoref, pszFilename );
329 : }
330 :
331 82 : 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 82 : if( fpBlob != NULL )
344 82 : VSIFCloseL( fpBlob );
345 :
346 82 : if( nGCPCount > 0 )
347 : {
348 22 : GDALDeinitGCPs( nGCPCount, pasGCPList );
349 22 : CPLFree( pasGCPList );
350 : }
351 :
352 82 : CPLFree( pszProjection );
353 82 : CPLFree( pszGCPProjection );
354 82 : CPLFree( pszPath );
355 82 : CSLDestroy( papszGeoref );
356 82 : CSLDestroy( papszAttrib );
357 82 : }
358 :
359 : /************************************************************************/
360 : /* SetVersion() */
361 : /************************************************************************/
362 :
363 82 : void HKVDataset::SetVersion(float version_number)
364 :
365 : {
366 : //update stored info
367 82 : MFF2version = version_number;
368 82 : }
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 50 : 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 50 : pszFilename = CPLFormFilename( pszFilenameIn, "attrib", NULL );
408 :
409 50 : fp = VSIFOpen( pszFilename, "wt" );
410 50 : 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 50 : fprintf( fp, "channel.enumeration = %d\n", nBands );
418 50 : fprintf( fp, "channel.interleave = { *pixel tile sequential }\n" );
419 50 : fprintf( fp, "extent.cols = %d\n", nXSize );
420 50 : fprintf( fp, "extent.rows = %d\n", nYSize );
421 :
422 50 : switch( eType )
423 : {
424 : case GDT_Byte:
425 : fprintf( fp, "pixel.encoding = "
426 20 : "{ *unsigned twos-complement ieee-754 }\n" );
427 20 : break;
428 :
429 : case GDT_UInt16:
430 : fprintf( fp, "pixel.encoding = "
431 6 : "{ *unsigned twos-complement ieee-754 }\n" );
432 6 : break;
433 :
434 : case GDT_CInt16:
435 : case GDT_Int16:
436 : fprintf( fp, "pixel.encoding = "
437 12 : "{ unsigned *twos-complement ieee-754 }\n" );
438 12 : break;
439 :
440 : case GDT_CFloat32:
441 : case GDT_Float32:
442 : fprintf( fp, "pixel.encoding = "
443 12 : "{ unsigned twos-complement *ieee-754 }\n" );
444 12 : break;
445 :
446 : default:
447 0 : CPLAssert( FALSE );
448 : }
449 :
450 50 : fprintf( fp, "pixel.size = %d\n", GDALGetDataTypeSize(eType) );
451 50 : if( GDALDataTypeIsComplex( eType ) )
452 12 : fprintf( fp, "pixel.field = { real *complex }\n" );
453 : else
454 38 : fprintf( fp, "pixel.field = { *real complex }\n" );
455 :
456 : #ifdef CPL_MSB
457 : fprintf( fp, "pixel.order = { lsbf *msbf }\n" );
458 : #else
459 50 : fprintf( fp, "pixel.order = { *lsbf msbf }\n" );
460 : #endif
461 :
462 50 : if ( bNoDataSet )
463 0 : fprintf( fp, "pixel.no_data = %f\n", dfNoDataValue );
464 :
465 : /* version information- only create the new style */
466 50 : fprintf( fp, "version = 1.1");
467 :
468 :
469 50 : VSIFClose( fp );
470 50 : return CE_None;
471 : }
472 :
473 :
474 : /************************************************************************/
475 : /* GetProjectionRef() */
476 : /************************************************************************/
477 :
478 40 : const char *HKVDataset::GetProjectionRef()
479 :
480 : {
481 40 : return( pszProjection );
482 : }
483 :
484 : /************************************************************************/
485 : /* GetGeoTransform() */
486 : /************************************************************************/
487 :
488 50 : CPLErr HKVDataset::GetGeoTransform( double * padfTransform )
489 :
490 : {
491 50 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
492 50 : return( CE_None );
493 : }
494 :
495 : /************************************************************************/
496 : /* SetGeoTransform() */
497 : /************************************************************************/
498 :
499 50 : 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 50 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
512 :
513 : /* Clear previous gcps */
514 50 : if( nGCPCount > 0 )
515 : {
516 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
517 0 : CPLFree( pasGCPList );
518 : }
519 50 : nGCPCount = 0;
520 50 : pasGCPList = NULL;
521 :
522 : /* Return if the identity transform is set */
523 50 : 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 50 : OGRSpatialReference oUTM;
533 50 : OGRSpatialReference oLL;
534 50 : OGRCoordinateTransformation *poTransform = NULL;
535 50 : int bSuccess=TRUE;
536 : char *pszPtemp;
537 : char *pszGCPtemp;
538 :
539 : /* Projection parameter checking will have been done */
540 : /* in SetProjection. */
541 50 : 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 50 : else if ((( CSLFetchNameValue( papszGeoref, "projection.name" ) != NULL ) &&
559 : ( !EQUAL(CSLFetchNameValue( papszGeoref, "projection.name" ),"LL" ))) ||
560 : (CSLFetchNameValue( papszGeoref, "projection.name" ) == NULL ))
561 : {
562 30 : return CE_Failure;
563 : }
564 :
565 20 : nGCPCount = 0;
566 20 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5);
567 :
568 : /* -------------------------------------------------------------------- */
569 : /* top left */
570 : /* -------------------------------------------------------------------- */
571 20 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
572 20 : CPLFree( pasGCPList[nGCPCount].pszId );
573 20 : pasGCPList[nGCPCount].pszId = CPLStrdup( "top_left" );
574 :
575 20 : if (MFF2version > 1.0)
576 : {
577 20 : temp_lat = padfTransform[3];
578 20 : temp_long = padfTransform[0];
579 20 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
580 20 : 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 20 : pasGCPList[nGCPCount].dfGCPX = temp_long;
590 20 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
591 20 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
592 20 : nGCPCount++;
593 :
594 20 : if (poTransform != NULL)
595 : {
596 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
597 0 : bSuccess = FALSE;
598 : }
599 :
600 20 : if (bSuccess)
601 : {
602 20 : sprintf( szValue, "%.10f", temp_lat );
603 : papszGeoref = CSLSetNameValue( papszGeoref, "top_left.latitude",
604 20 : szValue );
605 :
606 20 : sprintf( szValue, "%.10f", temp_long );
607 : papszGeoref = CSLSetNameValue( papszGeoref, "top_left.longitude",
608 20 : szValue );
609 : }
610 :
611 : /* -------------------------------------------------------------------- */
612 : /* top_right */
613 : /* -------------------------------------------------------------------- */
614 20 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
615 20 : CPLFree( pasGCPList[nGCPCount].pszId );
616 20 : pasGCPList[nGCPCount].pszId = CPLStrdup( "top_right" );
617 :
618 20 : if (MFF2version > 1.0)
619 : {
620 20 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4];
621 20 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1];
622 20 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
623 20 : 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 20 : pasGCPList[nGCPCount].dfGCPX = temp_long;
633 20 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
634 20 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
635 20 : nGCPCount++;
636 :
637 20 : if (poTransform != NULL)
638 : {
639 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
640 0 : bSuccess = FALSE;
641 : }
642 :
643 20 : if (bSuccess)
644 : {
645 20 : sprintf( szValue, "%.10f", temp_lat );
646 : papszGeoref = CSLSetNameValue( papszGeoref, "top_right.latitude",
647 20 : szValue );
648 :
649 20 : sprintf( szValue, "%.10f", temp_long );
650 : papszGeoref = CSLSetNameValue( papszGeoref, "top_right.longitude",
651 20 : szValue );
652 : }
653 :
654 : /* -------------------------------------------------------------------- */
655 : /* bottom_left */
656 : /* -------------------------------------------------------------------- */
657 20 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
658 20 : CPLFree( pasGCPList[nGCPCount].pszId );
659 20 : pasGCPList[nGCPCount].pszId = CPLStrdup( "bottom_left" );
660 :
661 20 : if (MFF2version > 1.0)
662 : {
663 20 : temp_lat = padfTransform[3] + GetRasterYSize() * padfTransform[5];
664 20 : temp_long = padfTransform[0] + GetRasterYSize() * padfTransform[2];
665 20 : pasGCPList[nGCPCount].dfGCPPixel = 0.0;
666 20 : 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 20 : pasGCPList[nGCPCount].dfGCPX = temp_long;
676 20 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
677 20 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
678 20 : nGCPCount++;
679 :
680 20 : if (poTransform != NULL)
681 : {
682 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
683 0 : bSuccess = FALSE;
684 : }
685 :
686 20 : if (bSuccess)
687 : {
688 20 : sprintf( szValue, "%.10f", temp_lat );
689 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_left.latitude",
690 20 : szValue );
691 :
692 20 : sprintf( szValue, "%.10f", temp_long );
693 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_left.longitude",
694 20 : szValue );
695 : }
696 :
697 : /* -------------------------------------------------------------------- */
698 : /* bottom_right */
699 : /* -------------------------------------------------------------------- */
700 20 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
701 20 : CPLFree( pasGCPList[nGCPCount].pszId );
702 20 : pasGCPList[nGCPCount].pszId = CPLStrdup( "bottom_right" );
703 :
704 20 : if (MFF2version > 1.0)
705 : {
706 40 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] +
707 40 : GetRasterYSize() * padfTransform[5];
708 40 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] +
709 40 : GetRasterYSize() * padfTransform[2];
710 20 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize();
711 20 : 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 20 : pasGCPList[nGCPCount].dfGCPX = temp_long;
724 20 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
725 20 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
726 20 : nGCPCount++;
727 :
728 20 : if (poTransform != NULL)
729 : {
730 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
731 0 : bSuccess = FALSE;
732 : }
733 :
734 20 : if (bSuccess)
735 : {
736 20 : sprintf( szValue, "%.10f", temp_lat );
737 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_right.latitude",
738 20 : szValue );
739 :
740 20 : sprintf( szValue, "%.10f", temp_long );
741 : papszGeoref = CSLSetNameValue( papszGeoref, "bottom_right.longitude",
742 20 : szValue );
743 : }
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Center */
747 : /* -------------------------------------------------------------------- */
748 20 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
749 20 : CPLFree( pasGCPList[nGCPCount].pszId );
750 20 : pasGCPList[nGCPCount].pszId = CPLStrdup( "centre" );
751 :
752 20 : if (MFF2version > 1.0)
753 : {
754 40 : temp_lat = padfTransform[3] + GetRasterXSize() * padfTransform[4] * 0.5 +
755 40 : GetRasterYSize() * padfTransform[5] * 0.5;
756 40 : temp_long = padfTransform[0] + GetRasterXSize() * padfTransform[1] * 0.5 +
757 40 : GetRasterYSize() * padfTransform[2] * 0.5;
758 20 : pasGCPList[nGCPCount].dfGCPPixel = GetRasterXSize()/2.0;
759 20 : 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 20 : pasGCPList[nGCPCount].dfGCPX = temp_long;
771 20 : pasGCPList[nGCPCount].dfGCPY = temp_lat;
772 20 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
773 20 : nGCPCount++;
774 :
775 20 : if (poTransform != NULL)
776 : {
777 0 : if( !bSuccess || !poTransform->Transform( 1, &temp_long, &temp_lat ) )
778 0 : bSuccess = FALSE;
779 : }
780 :
781 20 : if (bSuccess)
782 : {
783 20 : sprintf( szValue, "%.10f", temp_lat );
784 : papszGeoref = CSLSetNameValue( papszGeoref, "centre.latitude",
785 20 : szValue );
786 :
787 20 : sprintf( szValue, "%.10f", temp_long );
788 : papszGeoref = CSLSetNameValue( papszGeoref, "centre.longitude",
789 20 : szValue );
790 : }
791 :
792 20 : 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 20 : if (poTransform != NULL)
799 0 : delete poTransform;
800 :
801 :
802 :
803 20 : bGeorefChanged = TRUE;
804 :
805 20 : return( CE_None );
806 : }
807 :
808 20 : CPLErr HKVDataset::SetGCPProjection( const char *pszNewProjection )
809 : {
810 :
811 20 : CPLFree( pszGCPProjection );
812 20 : this->pszGCPProjection = CPLStrdup(pszNewProjection);
813 :
814 20 : return CE_None;
815 : }
816 :
817 : /************************************************************************/
818 : /* SetProjection() */
819 : /* */
820 : /* We provide very limited support for setting the projection. */
821 : /************************************************************************/
822 :
823 50 : CPLErr HKVDataset::SetProjection( const char * pszNewProjection )
824 :
825 : {
826 : HKVSpheroidList *hkvEllipsoids;
827 : double eq_radius, inv_flattening;
828 50 : OGRErr ogrerrorEq=OGRERR_NONE;
829 50 : OGRErr ogrerrorInvf=OGRERR_NONE;
830 50 : OGRErr ogrerrorOl=OGRERR_NONE;
831 :
832 50 : 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 50 : 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 50 : else if (EQUAL(pszNewProjection,""))
851 : {
852 0 : CPLFree( pszProjection );
853 0 : pszProjection = (char *) CPLStrdup(pszNewProjection);
854 :
855 0 : return CE_None;
856 : }
857 50 : CPLFree( pszProjection );
858 50 : pszProjection = (char *) CPLStrdup(pszNewProjection);
859 :
860 :
861 50 : OGRSpatialReference oSRS(pszNewProjection);
862 :
863 50 : 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 50 : else if ((oSRS.GetAttrValue("PROJECTION") == NULL) && (oSRS.IsGeographic()))
875 : {
876 50 : 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 50 : eq_radius = oSRS.GetSemiMajor(&ogrerrorEq);
885 50 : inv_flattening = oSRS.GetInvFlattening(&ogrerrorInvf);
886 100 : if ((ogrerrorEq == OGRERR_NONE) && (ogrerrorInvf == OGRERR_NONE))
887 : {
888 50 : hkvEllipsoids = new HKVSpheroidList;
889 100 : spheroid_name = hkvEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(eq_radius,inv_flattening);
890 50 : if (spheroid_name != NULL)
891 : {
892 : papszGeoref = CSLSetNameValue( papszGeoref, "spheroid.name",
893 50 : spheroid_name );
894 : }
895 50 : CPLFree(spheroid_name);
896 50 : 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 50 : bGeorefChanged = TRUE;
915 50 : 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 160 : 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 160 : sprintf( szFieldName, "%s.latitude", pszBase );
963 160 : if( CSLFetchNameValue(papszGeoref, szFieldName) == NULL )
964 150 : return;
965 : else
966 10 : dfLat = atof(CSLFetchNameValue(papszGeoref, szFieldName));
967 :
968 10 : sprintf( szFieldName, "%s.longitude", pszBase );
969 10 : if( CSLFetchNameValue(papszGeoref, szFieldName) == NULL )
970 0 : return;
971 : else
972 10 : dfLong = atof(CSLFetchNameValue(papszGeoref, szFieldName));
973 :
974 : /* -------------------------------------------------------------------- */
975 : /* Add the gcp to the internal list. */
976 : /* -------------------------------------------------------------------- */
977 10 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
978 :
979 10 : CPLFree( pasGCPList[nGCPCount].pszId );
980 :
981 10 : pasGCPList[nGCPCount].pszId = CPLStrdup( pszBase );
982 :
983 10 : pasGCPList[nGCPCount].dfGCPX = dfLong;
984 10 : pasGCPList[nGCPCount].dfGCPY = dfLat;
985 10 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
986 :
987 10 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
988 10 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
989 :
990 10 : nGCPCount++;
991 : }
992 :
993 : /************************************************************************/
994 : /* ProcessGeoref() */
995 : /************************************************************************/
996 :
997 32 : void HKVDataset::ProcessGeoref( const char * pszFilename )
998 :
999 : {
1000 : int i;
1001 32 : HKVSpheroidList *hkvEllipsoids = NULL;
1002 :
1003 : /* -------------------------------------------------------------------- */
1004 : /* Load the georef file, and boil white space away from around */
1005 : /* the equal sign. */
1006 : /* -------------------------------------------------------------------- */
1007 32 : CSLDestroy( papszGeoref );
1008 32 : papszGeoref = CSLLoad( pszFilename );
1009 32 : if( papszGeoref == NULL )
1010 0 : return;
1011 :
1012 32 : hkvEllipsoids = new HKVSpheroidList;
1013 :
1014 118 : for( i = 0; papszGeoref[i] != NULL; i++ )
1015 : {
1016 86 : int bAfterEqual = FALSE;
1017 : int iSrc, iDst;
1018 86 : char *pszLine = papszGeoref[i];
1019 :
1020 2066 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
1021 : {
1022 1980 : if( bAfterEqual || pszLine[iSrc] != ' ' )
1023 : {
1024 1980 : pszLine[iDst++] = pszLine[iSrc];
1025 : }
1026 :
1027 1980 : if( iDst > 0 && pszLine[iDst-1] == '=' )
1028 86 : bAfterEqual = FALSE;
1029 : }
1030 86 : pszLine[iDst] = '\0';
1031 : }
1032 :
1033 : /* -------------------------------------------------------------------- */
1034 : /* Try to get GCPs, in lat/longs . */
1035 : /* -------------------------------------------------------------------- */
1036 32 : nGCPCount = 0;
1037 32 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5);
1038 :
1039 32 : if (MFF2version > 1.0)
1040 : {
1041 : ProcessGeorefGCP( papszGeoref, "top_left",
1042 32 : 0, 0 );
1043 : ProcessGeorefGCP( papszGeoref, "top_right",
1044 32 : GetRasterXSize(), 0 );
1045 : ProcessGeorefGCP( papszGeoref, "bottom_left",
1046 32 : 0, GetRasterYSize() );
1047 : ProcessGeorefGCP( papszGeoref, "bottom_right",
1048 32 : GetRasterXSize(), GetRasterYSize() );
1049 : ProcessGeorefGCP( papszGeoref, "centre",
1050 32 : 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 32 : if (nGCPCount == 0)
1067 : {
1068 30 : CPLFree(pasGCPList);
1069 30 : 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 32 : "projection.name");
1080 : pszOriginLong = CSLFetchNameValue(papszGeoref,
1081 32 : "projection.origin_longitude");
1082 : pszSpheroidName = CSLFetchNameValue(papszGeoref,
1083 32 : "spheroid.name");
1084 :
1085 :
1086 32 : if ((pszSpheroidName != NULL) && (hkvEllipsoids->SpheroidInList(pszSpheroidName)))
1087 : {
1088 32 : eq_radius=hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
1089 32 : 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 34 : if( (pszProjName != NULL) && EQUAL(pszProjName,"utm") && (nGCPCount == 5) )
1099 : {
1100 : /*int nZone = (int)((atof(pszOriginLong)+184.5) / 6.0); */
1101 : int nZone;
1102 :
1103 2 : 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 2 : nZone = 31 + (int) floor(atof(pszOriginLong)/6.0);
1112 :
1113 2 : OGRSpatialReference oUTM;
1114 2 : OGRSpatialReference oLL;
1115 2 : OGRCoordinateTransformation *poTransform = NULL;
1116 : double dfUtmX[5], dfUtmY[5];
1117 : int gcp_index;
1118 :
1119 2 : int bSuccess = TRUE;
1120 :
1121 2 : if( pasGCPList[4].dfGCPY < 0 )
1122 0 : oUTM.SetUTM( nZone, 0 );
1123 : else
1124 2 : oUTM.SetUTM( nZone, 1 );
1125 :
1126 2 : if (pszOriginLong != NULL)
1127 : {
1128 2 : oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN,atof(pszOriginLong));
1129 2 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));
1130 : }
1131 :
1132 2 : 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 2 : if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
1141 : {
1142 : oUTM.SetGeogCS( "unknown","unknown",pszSpheroidName,
1143 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1144 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1145 2 : );
1146 : oLL.SetGeogCS( "unknown","unknown",pszSpheroidName,
1147 : hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
1148 : hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
1149 2 : );
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 2 : poTransform = OGRCreateCoordinateTransformation( &oLL, &oUTM );
1160 2 : if( poTransform == NULL )
1161 : {
1162 0 : CPLErrorReset();
1163 0 : bSuccess = FALSE;
1164 : }
1165 :
1166 12 : for(gcp_index=0;gcp_index<5;gcp_index++)
1167 : {
1168 10 : dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
1169 10 : dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
1170 :
1171 10 : if( bSuccess && !poTransform->Transform( 1, &(dfUtmX[gcp_index]), &(dfUtmY[gcp_index]) ) )
1172 0 : bSuccess = FALSE;
1173 :
1174 : }
1175 :
1176 2 : if( bSuccess )
1177 : {
1178 2 : int transform_ok = FALSE;
1179 :
1180 : /* update GCPS to proper projection */
1181 12 : for(gcp_index=0;gcp_index<5;gcp_index++)
1182 : {
1183 10 : pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
1184 10 : pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
1185 : }
1186 :
1187 2 : CPLFree( pszGCPProjection );
1188 2 : pszGCPProjection = NULL;
1189 2 : oUTM.exportToWkt( &pszGCPProjection );
1190 :
1191 2 : transform_ok = GDALGCPsToGeoTransform(5,pasGCPList,adfGeoTransform,0);
1192 :
1193 2 : CPLFree( pszProjection );
1194 2 : pszProjection = NULL;
1195 2 : 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 2 : oUTM.exportToWkt( &pszProjection );
1208 :
1209 : }
1210 :
1211 2 : if( poTransform != NULL )
1212 2 : delete poTransform;
1213 : }
1214 30 : 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 32 : delete hkvEllipsoids;
1272 : }
1273 :
1274 : /************************************************************************/
1275 : /* Open() */
1276 : /************************************************************************/
1277 :
1278 22746 : GDALDataset *HKVDataset::Open( GDALOpenInfo * poOpenInfo )
1279 :
1280 : {
1281 22746 : int i, bNoDataSet = FALSE;
1282 22746 : 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 22746 : if( !poOpenInfo->bIsDirectory )
1292 22662 : return NULL;
1293 :
1294 84 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "image_data", NULL);
1295 84 : if( VSIStat(pszFilename,&sStat) != 0 )
1296 2 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "blob", NULL );
1297 84 : if( VSIStat(pszFilename,&sStat) != 0 )
1298 2 : return NULL;
1299 :
1300 82 : pszFilename = CPLFormFilename(poOpenInfo->pszFilename, "attrib", NULL );
1301 82 : 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 82 : papszAttrib = CSLLoad( pszFilename );
1309 82 : if( papszAttrib == NULL )
1310 0 : return NULL;
1311 :
1312 820 : for( i = 0; papszAttrib[i] != NULL; i++ )
1313 : {
1314 738 : int bAfterEqual = FALSE;
1315 : int iSrc, iDst;
1316 738 : char *pszLine = papszAttrib[i];
1317 :
1318 20834 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
1319 : {
1320 20096 : if( bAfterEqual || pszLine[iSrc] != ' ' )
1321 : {
1322 17472 : pszLine[iDst++] = pszLine[iSrc];
1323 : }
1324 :
1325 20096 : if( iDst > 0 && pszLine[iDst-1] == '=' )
1326 1476 : bAfterEqual = FALSE;
1327 : }
1328 738 : pszLine[iDst] = '\0';
1329 : }
1330 :
1331 : /* -------------------------------------------------------------------- */
1332 : /* Create a corresponding GDALDataset. */
1333 : /* -------------------------------------------------------------------- */
1334 : HKVDataset *poDS;
1335 :
1336 82 : poDS = new HKVDataset();
1337 :
1338 82 : poDS->pszPath = CPLStrdup( poOpenInfo->pszFilename );
1339 82 : poDS->papszAttrib = papszAttrib;
1340 :
1341 82 : poDS->eAccess = poOpenInfo->eAccess;
1342 :
1343 : /* -------------------------------------------------------------------- */
1344 : /* Set some dataset wide information. */
1345 : /* -------------------------------------------------------------------- */
1346 : int bNative, bComplex;
1347 82 : int nRawBands = 0;
1348 :
1349 164 : if( CSLFetchNameValue( papszAttrib, "extent.cols" ) == NULL
1350 : || CSLFetchNameValue( papszAttrib, "extent.rows" ) == NULL )
1351 0 : return NULL;
1352 :
1353 82 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib,"extent.cols"));
1354 82 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib,"extent.rows"));
1355 :
1356 82 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
1357 : {
1358 0 : delete poDS;
1359 0 : return NULL;
1360 : }
1361 :
1362 82 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.order");
1363 82 : if( pszValue == NULL )
1364 0 : bNative = TRUE;
1365 : else
1366 : {
1367 : #ifdef CPL_MSB
1368 : bNative = (strstr(pszValue,"*msbf") != NULL);
1369 : #else
1370 82 : bNative = (strstr(pszValue,"*lsbf") != NULL);
1371 : #endif
1372 : }
1373 :
1374 82 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.no_data");
1375 82 : if( pszValue != NULL )
1376 : {
1377 0 : bNoDataSet = TRUE;
1378 0 : dfNoDataValue = atof(pszValue);
1379 : }
1380 :
1381 82 : pszValue = CSLFetchNameValue(papszAttrib,"channel.enumeration");
1382 82 : if( pszValue != NULL )
1383 82 : nRawBands = atoi(pszValue);
1384 : else
1385 0 : nRawBands = 1;
1386 :
1387 82 : if (!GDALCheckBandCount(nRawBands, TRUE))
1388 : {
1389 0 : delete poDS;
1390 0 : return NULL;
1391 : }
1392 :
1393 82 : pszValue = CSLFetchNameValue(papszAttrib,"pixel.field");
1394 102 : if( pszValue != NULL && strstr(pszValue,"*complex") != NULL )
1395 20 : bComplex = TRUE;
1396 : else
1397 62 : 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 82 : if (CSLFetchNameValue( papszAttrib, "version" ) != NULL)
1403 : poDS->SetVersion((float)
1404 82 : 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 82 : int nSize = 1;
1413 : int nPseudoBands;
1414 : GDALDataType eType;
1415 :
1416 82 : pszEncoding = CSLFetchNameValue(papszAttrib,"pixel.encoding");
1417 82 : if( pszEncoding == NULL )
1418 0 : pszEncoding = "{ *unsigned }";
1419 :
1420 82 : if( CSLFetchNameValue(papszAttrib,"pixel.size") != NULL )
1421 82 : nSize = atoi(CSLFetchNameValue(papszAttrib,"pixel.size"))/8;
1422 :
1423 82 : if( bComplex )
1424 20 : nPseudoBands = 2;
1425 : else
1426 62 : nPseudoBands = 1;
1427 :
1428 82 : if( nSize == 1 )
1429 32 : eType = GDT_Byte;
1430 60 : else if( nSize == 2 && strstr(pszEncoding,"*unsigned") != NULL )
1431 10 : eType = GDT_UInt16;
1432 50 : else if( nSize == 4 && bComplex )
1433 10 : eType = GDT_CInt16;
1434 30 : else if( nSize == 2 )
1435 10 : eType = GDT_Int16;
1436 20 : else if( nSize == 4 && strstr(pszEncoding,"*unsigned") != NULL )
1437 0 : eType = GDT_UInt32;
1438 20 : else if( nSize == 8 && strstr(pszEncoding,"*two") != NULL && bComplex )
1439 0 : eType = GDT_CInt32;
1440 20 : else if( nSize == 4 && strstr(pszEncoding,"*two") != NULL )
1441 0 : eType = GDT_Int32;
1442 30 : else if( nSize == 8 && bComplex )
1443 10 : eType = GDT_CFloat32;
1444 10 : else if( nSize == 4 )
1445 10 : 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 82 : pszFilename = CPLFormFilename(poDS->pszPath, "image_data", NULL );
1464 82 : if( VSIStat(pszFilename,&sStat) != 0 )
1465 0 : pszFilename = CPLFormFilename(poDS->pszPath, "blob", NULL );
1466 82 : if( poOpenInfo->eAccess == GA_ReadOnly )
1467 : {
1468 32 : poDS->fpBlob = VSIFOpenL( pszFilename, "rb" );
1469 32 : 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 50 : poDS->fpBlob = VSIFOpenL( pszFilename, "rb+" );
1481 50 : 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 82 : char *pszOvrFilename = (char *) CPLMalloc(strlen(pszFilename)+5);
1495 :
1496 82 : sprintf( pszOvrFilename, "%s_ovr", pszFilename );
1497 :
1498 : /* -------------------------------------------------------------------- */
1499 : /* Define the bands. */
1500 : /* -------------------------------------------------------------------- */
1501 : int nPixelOffset, nLineOffset, nOffset;
1502 :
1503 82 : nPixelOffset = nRawBands * nSize;
1504 82 : nLineOffset = nPixelOffset * poDS->GetRasterXSize();
1505 82 : nOffset = 0;
1506 :
1507 264 : 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 182 : eType, bNative );
1515 182 : poDS->SetBand( poDS->GetRasterCount()+1, poBand );
1516 182 : nOffset += GDALGetDataTypeSize( eType ) / 8;
1517 :
1518 182 : if( bNoDataSet )
1519 0 : poBand->SetNoDataValue( dfNoDataValue );
1520 : }
1521 :
1522 82 : poDS->eRasterType = eType;
1523 :
1524 : /* -------------------------------------------------------------------- */
1525 : /* Process the georef file if there is one. */
1526 : /* -------------------------------------------------------------------- */
1527 82 : pszFilename = CPLFormFilename(poDS->pszPath, "georef", NULL );
1528 82 : if( VSIStat(pszFilename,&sStat) == 0 )
1529 32 : poDS->ProcessGeoref(pszFilename);
1530 :
1531 : /* -------------------------------------------------------------------- */
1532 : /* Initialize any PAM information. */
1533 : /* -------------------------------------------------------------------- */
1534 82 : poDS->SetDescription( pszOvrFilename );
1535 82 : poDS->TryLoadXML();
1536 :
1537 : /* -------------------------------------------------------------------- */
1538 : /* Handle overviews. */
1539 : /* -------------------------------------------------------------------- */
1540 82 : poDS->oOvManager.Initialize( poDS, pszOvrFilename, NULL, TRUE );
1541 :
1542 82 : CPLFree( pszOvrFilename );
1543 :
1544 82 : return( poDS );
1545 : }
1546 :
1547 : /************************************************************************/
1548 : /* Create() */
1549 : /************************************************************************/
1550 :
1551 86 : 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 86 : if (nBands <= 0)
1561 : {
1562 : CPLError( CE_Failure, CPLE_NotSupported,
1563 2 : "HKV driver does not support %d bands.\n", nBands);
1564 2 : return NULL;
1565 : }
1566 :
1567 84 : 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 30 : GDALGetDataTypeName(eType) );
1576 :
1577 30 : 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 54 : if( strlen(CPLGetPath(pszFilenameIn)) == 0 )
1588 0 : pszBaseDir = CPLStrdup(".");
1589 : else
1590 54 : pszBaseDir = CPLStrdup(CPLGetPath(pszFilenameIn));
1591 :
1592 54 : 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 4 : pszBaseDir);
1598 4 : CPLFree( pszBaseDir );
1599 4 : return NULL;
1600 : }
1601 :
1602 50 : CPLFree( pszBaseDir );
1603 50 : pszBaseDir = NULL;
1604 :
1605 50 : 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 50 : nBands, eType, FALSE, 0.0 );
1620 :
1621 50 : 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 50 : pszFilename = CPLFormFilename( pszFilenameIn, "image_data", NULL );
1632 50 : fp = VSIFOpen( pszFilename, "wb" );
1633 50 : 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 50 : VSIFWrite( (void*)"", 1, 1, fp );
1641 50 : VSIFClose( fp );
1642 :
1643 : /* -------------------------------------------------------------------- */
1644 : /* Open the dataset normally. */
1645 : /* -------------------------------------------------------------------- */
1646 50 : 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 36 : 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 36 : int nBands = poSrcDS->GetRasterCount();
1724 36 : if (nBands == 0)
1725 : {
1726 : CPLError( CE_Failure, CPLE_NotSupported,
1727 2 : "HKV driver does not support source dataset with zero band.\n");
1728 2 : return NULL;
1729 : }
1730 :
1731 34 : eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1732 :
1733 34 : 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 54 : 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 34 : eType, papszOptions );
1749 :
1750 : /* Check that Create worked- return Null if it didn't */
1751 34 : if (poDS == NULL)
1752 14 : return NULL;
1753 :
1754 : /* -------------------------------------------------------------------- */
1755 : /* Copy the image data. */
1756 : /* -------------------------------------------------------------------- */
1757 20 : int nXSize = poDS->GetRasterXSize();
1758 20 : int nYSize = poDS->GetRasterYSize();
1759 : int nBlockXSize, nBlockYSize, nBlockTotal, nBlocksDone;
1760 :
1761 20 : poDS->GetRasterBand(1)->GetBlockSize( &nBlockXSize, &nBlockYSize );
1762 :
1763 : nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize)
1764 : * ((nYSize + nBlockYSize - 1) / nBlockYSize)
1765 20 : * poSrcDS->GetRasterCount();
1766 :
1767 20 : nBlocksDone = 0;
1768 60 : for( iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++ )
1769 : {
1770 40 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand+1 );
1771 40 : GDALRasterBand *poDstBand = poDS->GetRasterBand( iBand+1 );
1772 : int iYOffset, iXOffset;
1773 : void *pData;
1774 : CPLErr eErr;
1775 : int pbSuccess;
1776 40 : double dfSrcNoDataValue =0.0;
1777 :
1778 : /* Get nodata value, if relevant */
1779 40 : dfSrcNoDataValue = poSrcBand->GetNoDataValue( &pbSuccess );
1780 40 : if ( pbSuccess )
1781 0 : poDS->SetNoDataValue( dfSrcNoDataValue );
1782 :
1783 : pData = CPLMalloc(nBlockXSize * nBlockYSize
1784 40 : * GDALGetDataTypeSize(eType) / 8);
1785 :
1786 440 : for( iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize )
1787 : {
1788 800 : for( iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize )
1789 : {
1790 : int nTBXSize, nTBYSize;
1791 :
1792 400 : 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 400 : nTBXSize = MIN(nBlockXSize,nXSize-iXOffset);
1806 400 : nTBYSize = MIN(nBlockYSize,nYSize-iYOffset);
1807 :
1808 : eErr = poSrcBand->RasterIO( GF_Read,
1809 : iXOffset, iYOffset,
1810 : nTBXSize, nTBYSize,
1811 : pData, nTBXSize, nTBYSize,
1812 400 : eType, 0, 0 );
1813 400 : 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 400 : eType, 0, 0 );
1823 :
1824 400 : if( eErr != CE_None )
1825 : {
1826 0 : return NULL;
1827 : }
1828 : }
1829 : }
1830 :
1831 40 : 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 20 : double *tempGeoTransform=NULL;
1840 :
1841 20 : tempGeoTransform = (double *) CPLMalloc(6*sizeof(double));
1842 :
1843 40 : if (( poSrcDS->GetGeoTransform( tempGeoTransform ) == CE_None)
1844 20 : && (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 20 : poDS->SetGCPProjection(poSrcDS->GetProjectionRef());
1850 20 : poDS->SetProjection(poSrcDS->GetProjectionRef());
1851 20 : poDS->SetGeoTransform(tempGeoTransform);
1852 :
1853 20 : 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 60 : for( iBand = 0; iBand < poDS->GetRasterCount(); iBand++ )
1866 : {
1867 40 : RawRasterBand *poDstBand = (RawRasterBand *) poDS->GetRasterBand( iBand+1 );
1868 40 : poDstBand->FlushCache();
1869 : }
1870 :
1871 :
1872 20 : 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 20 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1885 :
1886 20 : return poDS;
1887 : }
1888 :
1889 :
1890 : /************************************************************************/
1891 : /* GDALRegister_HKV() */
1892 : /************************************************************************/
1893 :
1894 1135 : void GDALRegister_HKV()
1895 :
1896 : {
1897 : GDALDriver *poDriver;
1898 :
1899 1135 : if( GDALGetDriverByName( "MFF2" ) == NULL )
1900 : {
1901 1093 : poDriver = new GDALDriver();
1902 :
1903 1093 : poDriver->SetDescription( "MFF2" );
1904 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1905 1093 : "Vexcel MFF2 (HKV) Raster" );
1906 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1907 1093 : "frmt_mff2.html" );
1908 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1909 1093 : "Byte Int16 UInt16 Int32 UInt32 CInt16 CInt32 Float32 Float64 CFloat32 CFloat64" );
1910 :
1911 1093 : poDriver->pfnOpen = HKVDataset::Open;
1912 1093 : poDriver->pfnCreate = HKVDataset::Create;
1913 1093 : poDriver->pfnDelete = HKVDataset::Delete;
1914 1093 : poDriver->pfnCreateCopy = HKVDataset::CreateCopy;
1915 :
1916 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
1917 : }
1918 1135 : }
|
