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