Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: ISIS Version 2 Driver
4 : * Purpose: Implementation of ISIS2Dataset
5 : * Author: Trent Hare (thare@usgs.gov),
6 : * Robert Soricone (rsoricone@usgs.gov)
7 : * Ludovic Mercier (ludovic.mercier@gmail.com)
8 : * Frank Warmerdam (warmerdam@pobox.com)
9 : *
10 : * NOTE: Original code authored by Trent and Robert and placed in the public
11 : * domain as per US government policy. I have (within my rights) appropriated
12 : * it and placed it under the following license. This is not intended to
13 : * diminish Trent and Roberts contribution.
14 : ******************************************************************************
15 : * Copyright (c) 2006, Frank Warmerdam <warmerdam@pobox.com>
16 : * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
17 : *
18 : * SPDX-License-Identifier: MIT
19 : ****************************************************************************/
20 :
21 : constexpr int NULL1 = 0;
22 : constexpr int NULL2 = -32768;
23 : constexpr double NULL3 = -3.4028226550889044521e+38;
24 :
25 : #include "cpl_string.h"
26 : #include "gdal_frmts.h"
27 : #include "nasakeywordhandler.h"
28 : #include "ogr_spatialref.h"
29 : #include "rawdataset.h"
30 : #include "pdsdrivercore.h"
31 :
32 : /************************************************************************/
33 : /* ==================================================================== */
34 : /* ISISDataset version2 */
35 : /* ==================================================================== */
36 : /************************************************************************/
37 :
38 : class ISIS2Dataset final : public RawDataset
39 : {
40 : VSILFILE *fpImage; // image data file.
41 : CPLString osExternalCube;
42 :
43 : NASAKeywordHandler oKeywords;
44 :
45 : int bGotTransform;
46 : double adfGeoTransform[6];
47 :
48 : OGRSpatialReference m_oSRS{};
49 :
50 : int parse_label(const char *file, char *keyword, char *value);
51 : int strstrip(char instr[], char outstr[], int position);
52 :
53 : CPLString oTempResult;
54 :
55 : static void CleanString(CPLString &osInput);
56 :
57 : const char *GetKeyword(const char *pszPath, const char *pszDefault = "");
58 : const char *GetKeywordSub(const char *pszPath, int iSubscript,
59 : const char *pszDefault = "");
60 :
61 : CPLErr Close() override;
62 :
63 : public:
64 : ISIS2Dataset();
65 : virtual ~ISIS2Dataset();
66 :
67 : virtual CPLErr GetGeoTransform(double *padfTransform) override;
68 : const OGRSpatialReference *GetSpatialRef() const override;
69 :
70 : virtual char **GetFileList() override;
71 :
72 : static GDALDataset *Open(GDALOpenInfo *);
73 : };
74 :
75 : /************************************************************************/
76 : /* ISIS2Dataset() */
77 : /************************************************************************/
78 :
79 2 : ISIS2Dataset::ISIS2Dataset() : fpImage(nullptr), bGotTransform(FALSE)
80 : {
81 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
82 2 : adfGeoTransform[0] = 0.0;
83 2 : adfGeoTransform[1] = 1.0;
84 2 : adfGeoTransform[2] = 0.0;
85 2 : adfGeoTransform[3] = 0.0;
86 2 : adfGeoTransform[4] = 0.0;
87 2 : adfGeoTransform[5] = 1.0;
88 2 : }
89 :
90 : /************************************************************************/
91 : /* ~ISIS2Dataset() */
92 : /************************************************************************/
93 :
94 4 : ISIS2Dataset::~ISIS2Dataset()
95 :
96 : {
97 2 : ISIS2Dataset::Close();
98 4 : }
99 :
100 : /************************************************************************/
101 : /* Close() */
102 : /************************************************************************/
103 :
104 4 : CPLErr ISIS2Dataset::Close()
105 : {
106 4 : CPLErr eErr = CE_None;
107 4 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
108 : {
109 2 : if (ISIS2Dataset::FlushCache(true) != CE_None)
110 0 : eErr = CE_Failure;
111 2 : if (fpImage != nullptr)
112 : {
113 2 : if (VSIFCloseL(fpImage) != 0)
114 0 : eErr = CE_Failure;
115 : }
116 2 : if (GDALPamDataset::Close() != CE_None)
117 0 : eErr = CE_Failure;
118 : }
119 4 : return eErr;
120 : }
121 :
122 : /************************************************************************/
123 : /* GetFileList() */
124 : /************************************************************************/
125 :
126 1 : char **ISIS2Dataset::GetFileList()
127 :
128 : {
129 1 : char **papszFileList = GDALPamDataset::GetFileList();
130 :
131 1 : if (!osExternalCube.empty())
132 0 : papszFileList = CSLAddString(papszFileList, osExternalCube);
133 :
134 1 : return papszFileList;
135 : }
136 :
137 : /************************************************************************/
138 : /* GetSpatialRef() */
139 : /************************************************************************/
140 :
141 1 : const OGRSpatialReference *ISIS2Dataset::GetSpatialRef() const
142 : {
143 1 : if (!m_oSRS.IsEmpty())
144 1 : return &m_oSRS;
145 0 : return GDALPamDataset::GetSpatialRef();
146 : }
147 :
148 : /************************************************************************/
149 : /* GetGeoTransform() */
150 : /************************************************************************/
151 :
152 1 : CPLErr ISIS2Dataset::GetGeoTransform(double *padfTransform)
153 :
154 : {
155 1 : if (bGotTransform)
156 : {
157 1 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
158 1 : return CE_None;
159 : }
160 :
161 0 : return GDALPamDataset::GetGeoTransform(padfTransform);
162 : }
163 :
164 : /************************************************************************/
165 : /* Open() */
166 : /************************************************************************/
167 :
168 2 : GDALDataset *ISIS2Dataset::Open(GDALOpenInfo *poOpenInfo)
169 : {
170 : /* -------------------------------------------------------------------- */
171 : /* Does this look like a CUBE or an IMAGE Primary Data Object? */
172 : /* -------------------------------------------------------------------- */
173 2 : if (!ISIS2DriverIdentify(poOpenInfo) || poOpenInfo->fpL == nullptr)
174 0 : return nullptr;
175 :
176 2 : VSILFILE *fpQube = poOpenInfo->fpL;
177 2 : poOpenInfo->fpL = nullptr;
178 :
179 4 : auto poDS = std::make_unique<ISIS2Dataset>();
180 :
181 2 : if (!poDS->oKeywords.Ingest(fpQube, 0))
182 : {
183 0 : VSIFCloseL(fpQube);
184 0 : return nullptr;
185 : }
186 :
187 2 : VSIFCloseL(fpQube);
188 :
189 : /* -------------------------------------------------------------------- */
190 : /* We assume the user is pointing to the label (i.e. .lab) file. */
191 : /* -------------------------------------------------------------------- */
192 : // QUBE can be inline or detached and point to an image name
193 : // ^QUBE = 76
194 : // ^QUBE = ("ui31s015.img",6441<BYTES>) - has another label on the image
195 : // ^QUBE = "ui31s015.img" - which implies no label or skip value
196 :
197 2 : const char *pszQube = poDS->GetKeyword("^QUBE");
198 2 : int nQube = 0;
199 2 : int bByteLocation = FALSE;
200 4 : CPLString osTargetFile = poOpenInfo->pszFilename;
201 :
202 2 : if (pszQube[0] == '"')
203 : {
204 0 : const CPLString osTPath = CPLGetPathSafe(poOpenInfo->pszFilename);
205 0 : CPLString osFilename = pszQube;
206 0 : poDS->CleanString(osFilename);
207 0 : osTargetFile = CPLFormCIFilenameSafe(osTPath, osFilename, nullptr);
208 0 : poDS->osExternalCube = osTargetFile;
209 : }
210 2 : else if (pszQube[0] == '(')
211 : {
212 0 : const CPLString osTPath = CPLGetPathSafe(poOpenInfo->pszFilename);
213 0 : CPLString osFilename = poDS->GetKeywordSub("^QUBE", 1, "");
214 0 : poDS->CleanString(osFilename);
215 0 : osTargetFile = CPLFormCIFilenameSafe(osTPath, osFilename, nullptr);
216 0 : poDS->osExternalCube = osTargetFile;
217 :
218 0 : nQube = atoi(poDS->GetKeywordSub("^QUBE", 2, "1"));
219 0 : if (strstr(poDS->GetKeywordSub("^QUBE", 2, "1"), "<BYTES>") != nullptr)
220 0 : bByteLocation = true;
221 : }
222 : else
223 : {
224 2 : nQube = atoi(pszQube);
225 2 : if (strstr(pszQube, "<BYTES>") != nullptr)
226 0 : bByteLocation = true;
227 : }
228 :
229 : /* -------------------------------------------------------------------- */
230 : /* Check if file an ISIS2 header file? Read a few lines of text */
231 : /* searching for something starting with nrows or ncols. */
232 : /* -------------------------------------------------------------------- */
233 :
234 : /* -------------------------------------------------------------------- */
235 : /* Checks to see if this is valid ISIS2 cube */
236 : /* SUFFIX_ITEM tag in .cub file should be (0,0,0); no side-planes */
237 : /* -------------------------------------------------------------------- */
238 2 : const int s_ix = atoi(poDS->GetKeywordSub("QUBE.SUFFIX_ITEMS", 1));
239 2 : const int s_iy = atoi(poDS->GetKeywordSub("QUBE.SUFFIX_ITEMS", 2));
240 2 : const int s_iz = atoi(poDS->GetKeywordSub("QUBE.SUFFIX_ITEMS", 3));
241 :
242 2 : if (s_ix != 0 || s_iy != 0 || s_iz != 0)
243 : {
244 0 : CPLError(CE_Failure, CPLE_OpenFailed,
245 : "*** ISIS 2 cube file has invalid SUFFIX_ITEMS parameters:\n"
246 : "*** gdal isis2 driver requires (0, 0, 0), thus no sideplanes "
247 : "or backplanes\n"
248 : "found: (%i, %i, %i)\n\n",
249 : s_ix, s_iy, s_iz);
250 0 : return nullptr;
251 : }
252 :
253 : /**************** end SUFFIX_ITEM check ***********************/
254 :
255 : /*********** Grab layout type (BSQ, BIP, BIL) ************/
256 : // AXIS_NAME = (SAMPLE,LINE,BAND)
257 : /***********************************************************/
258 :
259 2 : char szLayout[10] = "BSQ"; // default to band seq.
260 2 : const char *value = poDS->GetKeyword("QUBE.AXIS_NAME", "");
261 2 : if (EQUAL(value, "(SAMPLE,LINE,BAND)"))
262 2 : strcpy(szLayout, "BSQ");
263 0 : else if (EQUAL(value, "(BAND,LINE,SAMPLE)"))
264 0 : strcpy(szLayout, "BIP");
265 0 : else if (EQUAL(value, "(SAMPLE,BAND,LINE)") || EQUAL(value, ""))
266 0 : strcpy(szLayout, "BSQ");
267 : else
268 : {
269 0 : CPLError(CE_Failure, CPLE_OpenFailed,
270 : "%s layout not supported. Abort\n\n", value);
271 0 : return nullptr;
272 : }
273 :
274 : /*********** Grab samples lines band ************/
275 2 : const int nCols = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS", 1));
276 2 : const int nRows = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS", 2));
277 2 : const int nBands = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS", 3));
278 :
279 : /*********** Grab Qube record bytes **********/
280 2 : const int record_bytes = atoi(poDS->GetKeyword("RECORD_BYTES"));
281 2 : if (record_bytes < 0)
282 : {
283 0 : return nullptr;
284 : }
285 :
286 2 : GUIntBig nSkipBytes = 0;
287 2 : if (nQube > 0 && bByteLocation)
288 0 : nSkipBytes = (nQube - 1);
289 2 : else if (nQube > 0)
290 2 : nSkipBytes = static_cast<GUIntBig>(nQube - 1) * record_bytes;
291 : else
292 0 : nSkipBytes = 0;
293 :
294 : /*********** Grab samples lines band ************/
295 2 : char chByteOrder = 'M'; // default to MSB
296 4 : CPLString osCoreItemType = poDS->GetKeyword("QUBE.CORE_ITEM_TYPE");
297 2 : if ((EQUAL(osCoreItemType, "PC_INTEGER")) ||
298 4 : (EQUAL(osCoreItemType, "PC_UNSIGNED_INTEGER")) ||
299 2 : (EQUAL(osCoreItemType, "PC_REAL")))
300 : {
301 0 : chByteOrder = 'I';
302 : }
303 :
304 : /******** Grab format type - isis2 only supports 8,16,32 *******/
305 2 : GDALDataType eDataType = GDT_Byte;
306 2 : bool bNoDataSet = false;
307 2 : double dfNoData = 0.0;
308 :
309 2 : int itype = atoi(poDS->GetKeyword("QUBE.CORE_ITEM_BYTES", ""));
310 2 : switch (itype)
311 : {
312 0 : case 1:
313 0 : eDataType = GDT_Byte;
314 0 : dfNoData = NULL1;
315 0 : bNoDataSet = true;
316 0 : break;
317 0 : case 2:
318 0 : if (strstr(osCoreItemType, "UNSIGNED") != nullptr)
319 : {
320 0 : dfNoData = 0;
321 0 : eDataType = GDT_UInt16;
322 : }
323 : else
324 : {
325 0 : dfNoData = NULL2;
326 0 : eDataType = GDT_Int16;
327 : }
328 0 : bNoDataSet = true;
329 0 : break;
330 2 : case 4:
331 2 : eDataType = GDT_Float32;
332 2 : dfNoData = NULL3;
333 2 : bNoDataSet = true;
334 2 : break;
335 0 : case 8:
336 0 : eDataType = GDT_Float64;
337 0 : dfNoData = NULL3;
338 0 : bNoDataSet = true;
339 0 : break;
340 0 : default:
341 0 : CPLError(CE_Failure, CPLE_AppDefined,
342 : "Itype of %d is not supported in ISIS 2.", itype);
343 0 : return nullptr;
344 : }
345 :
346 : /*********** Grab Cellsize ************/
347 2 : double dfXDim = 1.0;
348 2 : double dfYDim = 1.0;
349 :
350 2 : value = poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.MAP_SCALE");
351 2 : if (strlen(value) > 0)
352 : {
353 : // Convert km to m
354 2 : dfXDim = static_cast<float>(CPLAtof(value) * 1000.0);
355 2 : dfYDim = static_cast<float>(CPLAtof(value) * 1000.0 * -1);
356 : }
357 :
358 : /*********** Grab LINE_PROJECTION_OFFSET ************/
359 2 : double dfULYMap = 0.5;
360 :
361 : value =
362 2 : poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.LINE_PROJECTION_OFFSET");
363 2 : if (strlen(value) > 0)
364 : {
365 2 : const double yulcenter = static_cast<float>(CPLAtof(value)) * dfYDim;
366 2 : dfULYMap = yulcenter - (dfYDim / 2);
367 : }
368 :
369 : /*********** Grab SAMPLE_PROJECTION_OFFSET ************/
370 2 : double dfULXMap = 0.5;
371 :
372 : value =
373 2 : poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.SAMPLE_PROJECTION_OFFSET");
374 2 : if (strlen(value) > 0)
375 : {
376 2 : const double xulcenter = static_cast<float>(CPLAtof(value)) * dfXDim;
377 2 : dfULXMap = xulcenter - (dfXDim / 2);
378 : }
379 :
380 : /*********** Grab TARGET_NAME ************/
381 : /**** This is the planets name i.e. MARS ***/
382 4 : const CPLString target_name = poDS->GetKeyword("QUBE.TARGET_NAME");
383 :
384 : /*********** Grab MAP_PROJECTION_TYPE ************/
385 : CPLString map_proj_name =
386 4 : poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.MAP_PROJECTION_TYPE");
387 2 : poDS->CleanString(map_proj_name);
388 :
389 : /*********** Grab SEMI-MAJOR ************/
390 : const double semi_major =
391 2 : CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.A_AXIS_RADIUS")) *
392 2 : 1000.0;
393 :
394 : /*********** Grab semi-minor ************/
395 : const double semi_minor =
396 2 : CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.C_AXIS_RADIUS")) *
397 2 : 1000.0;
398 :
399 : /*********** Grab CENTER_LAT ************/
400 : const double center_lat =
401 2 : CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.CENTER_LATITUDE"));
402 :
403 : /*********** Grab CENTER_LON ************/
404 : const double center_lon =
405 2 : CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.CENTER_LONGITUDE"));
406 :
407 : /*********** Grab 1st std parallel ************/
408 2 : const double first_std_parallel = CPLAtof(
409 : poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.FIRST_STANDARD_PARALLEL"));
410 :
411 : /*********** Grab 2nd std parallel ************/
412 2 : const double second_std_parallel = CPLAtof(
413 : poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.SECOND_STANDARD_PARALLEL"));
414 :
415 : /*** grab PROJECTION_LATITUDE_TYPE = "PLANETOCENTRIC" ****/
416 : // Need to further study how ocentric/ographic will effect the gdal library.
417 : // So far we will use this fact to define a sphere or ellipse for some
418 : // projections Frank - may need to talk this over
419 2 : bool bIsGeographic = true;
420 : value =
421 2 : poDS->GetKeyword("CUBE.IMAGE_MAP_PROJECTION.PROJECTION_LATITUDE_TYPE");
422 2 : if (EQUAL(value, "\"PLANETOCENTRIC\""))
423 0 : bIsGeographic = false;
424 :
425 2 : CPLDebug("ISIS2", "using projection %s", map_proj_name.c_str());
426 :
427 4 : OGRSpatialReference oSRS;
428 2 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
429 2 : bool bProjectionSet = true;
430 :
431 : // Set oSRS projection and parameters
432 2 : if ((EQUAL(map_proj_name, "EQUIRECTANGULAR_CYLINDRICAL")) ||
433 4 : (EQUAL(map_proj_name, "EQUIRECTANGULAR")) ||
434 2 : (EQUAL(map_proj_name, "SIMPLE_CYLINDRICAL")))
435 : {
436 2 : oSRS.OGRSpatialReference::SetEquirectangular2(0.0, center_lon,
437 : center_lat, 0, 0);
438 : }
439 0 : else if (EQUAL(map_proj_name, "ORTHOGRAPHIC"))
440 : {
441 0 : oSRS.OGRSpatialReference::SetOrthographic(center_lat, center_lon, 0, 0);
442 : }
443 0 : else if ((EQUAL(map_proj_name, "SINUSOIDAL")) ||
444 0 : (EQUAL(map_proj_name, "SINUSOIDAL_EQUAL-AREA")))
445 : {
446 0 : oSRS.OGRSpatialReference::SetSinusoidal(center_lon, 0, 0);
447 : }
448 0 : else if (EQUAL(map_proj_name, "MERCATOR"))
449 : {
450 0 : oSRS.OGRSpatialReference::SetMercator(center_lat, center_lon, 1, 0, 0);
451 : }
452 0 : else if (EQUAL(map_proj_name, "POLAR_STEREOGRAPHIC"))
453 : {
454 0 : oSRS.OGRSpatialReference::SetPS(center_lat, center_lon, 1, 0, 0);
455 : }
456 0 : else if (EQUAL(map_proj_name, "TRANSVERSE_MERCATOR"))
457 : {
458 0 : oSRS.OGRSpatialReference::SetTM(center_lat, center_lon, 1, 0, 0);
459 : }
460 0 : else if (EQUAL(map_proj_name, "LAMBERT_CONFORMAL_CONIC"))
461 : {
462 0 : oSRS.OGRSpatialReference::SetLCC(first_std_parallel,
463 : second_std_parallel, center_lat,
464 : center_lon, 0, 0);
465 : }
466 0 : else if (EQUAL(map_proj_name, ""))
467 : {
468 : /* no projection */
469 0 : bProjectionSet = false;
470 : }
471 : else
472 : {
473 0 : CPLDebug("ISIS2",
474 : "Dataset projection %s is not supported. Continuing...",
475 : map_proj_name.c_str());
476 0 : bProjectionSet = false;
477 : }
478 :
479 2 : if (bProjectionSet)
480 : {
481 : // Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
482 6 : const CPLString proj_target_name = map_proj_name + " " + target_name;
483 2 : oSRS.SetProjCS(proj_target_name); // set ProjCS keyword
484 :
485 : // The geographic/geocentric name will be the same basic name as the
486 : // body name 'GCS' = Geographic/Geocentric Coordinate System
487 4 : const CPLString geog_name = "GCS_" + target_name;
488 :
489 : // The datum and sphere names will be the same basic name aas the planet
490 4 : const CPLString datum_name = "D_" + target_name;
491 : // Might not be IAU defined so don't add.
492 4 : CPLString sphere_name = target_name; // + "_IAU_IAG");
493 :
494 : // calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
495 2 : double iflattening = 0.0;
496 2 : if ((semi_major - semi_minor) < 0.0000001)
497 2 : iflattening = 0;
498 : else
499 0 : iflattening = semi_major / (semi_major - semi_minor);
500 :
501 : // Set the body size but take into consideration which proj is being
502 : // used to help w/ proj4 compatibility The use of a Sphere, polar radius
503 : // or ellipse here is based on how ISIS does it internally
504 2 : if (((EQUAL(map_proj_name, "STEREOGRAPHIC") &&
505 4 : (fabs(center_lat) == 90))) ||
506 2 : (EQUAL(map_proj_name, "POLAR_STEREOGRAPHIC")))
507 : {
508 0 : if (bIsGeographic)
509 : {
510 : // Geograpraphic, so set an ellipse
511 0 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
512 : iflattening, "Reference_Meridian", 0.0);
513 : }
514 : else
515 : {
516 : // Geocentric, so force a sphere using the semi-minor axis. I
517 : // hope...
518 0 : sphere_name += "_polarRadius";
519 0 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_minor,
520 : 0.0, "Reference_Meridian", 0.0);
521 : }
522 : }
523 2 : else if ((EQUAL(map_proj_name, "SIMPLE_CYLINDRICAL")) ||
524 0 : (EQUAL(map_proj_name, "ORTHOGRAPHIC")) ||
525 0 : (EQUAL(map_proj_name, "STEREOGRAPHIC")) ||
526 2 : (EQUAL(map_proj_name, "SINUSOIDAL_EQUAL-AREA")) ||
527 0 : (EQUAL(map_proj_name, "SINUSOIDAL")))
528 : {
529 : // ISIS uses the spherical equation for these projections so force
530 : // a sphere.
531 2 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major, 0.0,
532 : "Reference_Meridian", 0.0);
533 : }
534 0 : else if ((EQUAL(map_proj_name, "EQUIRECTANGULAR_CYLINDRICAL")) ||
535 0 : (EQUAL(map_proj_name, "EQUIRECTANGULAR")))
536 : {
537 : // Calculate localRadius using ISIS3 simple elliptical method
538 : // not the more standard Radius of Curvature method
539 : // PI = 4 * atan(1);
540 0 : if (center_lon == 0)
541 : { // No need to calculate local radius
542 0 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
543 : 0.0, "Reference_Meridian", 0.0);
544 : }
545 : else
546 : {
547 0 : const double radLat = center_lat * M_PI / 180; // in radians
548 : const double localRadius =
549 0 : semi_major * semi_minor /
550 0 : sqrt(pow(semi_minor * cos(radLat), 2) +
551 0 : pow(semi_major * sin(radLat), 2));
552 0 : sphere_name += "_localRadius";
553 0 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, localRadius,
554 : 0.0, "Reference_Meridian", 0.0);
555 0 : CPLDebug("ISIS2", "local radius: %f", localRadius);
556 : }
557 : }
558 : else
559 : {
560 : // All other projections: Mercator, Transverse Mercator, Lambert
561 : // Conformal, etc. Geographic, so set an ellipse
562 0 : if (bIsGeographic)
563 : {
564 0 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
565 : iflattening, "Reference_Meridian", 0.0);
566 : }
567 : else
568 : {
569 : // Geocentric, so force a sphere. I hope...
570 0 : oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
571 : 0.0, "Reference_Meridian", 0.0);
572 : }
573 : }
574 :
575 : // translate back into a projection string.
576 2 : poDS->m_oSRS = std::move(oSRS);
577 : }
578 :
579 : /* END ISIS2 Label Read */
580 : /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
581 :
582 : /* -------------------------------------------------------------------- */
583 : /* Did we get the required keywords? If not we return with */
584 : /* this never having been considered to be a match. This isn't */
585 : /* an error! */
586 : /* -------------------------------------------------------------------- */
587 4 : if (!GDALCheckDatasetDimensions(nCols, nRows) ||
588 2 : !GDALCheckBandCount(nBands, false))
589 : {
590 0 : return nullptr;
591 : }
592 :
593 : /* -------------------------------------------------------------------- */
594 : /* Capture some information from the file that is of interest. */
595 : /* -------------------------------------------------------------------- */
596 2 : poDS->nRasterXSize = nCols;
597 2 : poDS->nRasterYSize = nRows;
598 :
599 : /* -------------------------------------------------------------------- */
600 : /* Open target binary file. */
601 : /* -------------------------------------------------------------------- */
602 :
603 2 : if (poOpenInfo->eAccess == GA_ReadOnly)
604 2 : poDS->fpImage = VSIFOpenL(osTargetFile, "rb");
605 : else
606 0 : poDS->fpImage = VSIFOpenL(osTargetFile, "r+b");
607 :
608 2 : if (poDS->fpImage == nullptr)
609 : {
610 0 : CPLError(CE_Failure, CPLE_OpenFailed,
611 : "Failed to open %s with write permission.\n%s",
612 0 : osTargetFile.c_str(), VSIStrerror(errno));
613 0 : return nullptr;
614 : }
615 :
616 2 : poDS->eAccess = poOpenInfo->eAccess;
617 :
618 : /* -------------------------------------------------------------------- */
619 : /* Compute the line offset. */
620 : /* -------------------------------------------------------------------- */
621 2 : int nItemSize = GDALGetDataTypeSizeBytes(eDataType);
622 : int nLineOffset, nPixelOffset;
623 : vsi_l_offset nBandOffset;
624 :
625 2 : if (EQUAL(szLayout, "BIP"))
626 : {
627 0 : nPixelOffset = nItemSize * nBands;
628 0 : if (nPixelOffset > INT_MAX / nBands)
629 : {
630 0 : return nullptr;
631 : }
632 0 : nLineOffset = nPixelOffset * nCols;
633 0 : nBandOffset = nItemSize;
634 : }
635 2 : else if (EQUAL(szLayout, "BSQ"))
636 : {
637 2 : nPixelOffset = nItemSize;
638 2 : if (nPixelOffset > INT_MAX / nCols)
639 : {
640 0 : return nullptr;
641 : }
642 2 : nLineOffset = nPixelOffset * nCols;
643 2 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
644 : }
645 : else /* assume BIL */
646 : {
647 0 : nPixelOffset = nItemSize;
648 0 : if (nPixelOffset > INT_MAX / nBands ||
649 0 : nPixelOffset * nBands > INT_MAX / nCols)
650 : {
651 0 : return nullptr;
652 : }
653 0 : nLineOffset = nItemSize * nBands * nCols;
654 0 : nBandOffset = static_cast<vsi_l_offset>(nItemSize) * nCols;
655 : }
656 :
657 : /* -------------------------------------------------------------------- */
658 : /* Create band information objects. */
659 : /* -------------------------------------------------------------------- */
660 4 : for (int i = 0; i < nBands; i++)
661 : {
662 : auto poBand = RawRasterBand::Create(
663 4 : poDS.get(), i + 1, poDS->fpImage, nSkipBytes + nBandOffset * i,
664 : nPixelOffset, nLineOffset, eDataType,
665 : chByteOrder == 'I' || chByteOrder == 'L'
666 2 : ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
667 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN,
668 2 : RawRasterBand::OwnFP::NO);
669 2 : if (!poBand)
670 0 : return nullptr;
671 :
672 2 : if (bNoDataSet)
673 2 : poBand->SetNoDataValue(dfNoData);
674 :
675 : // Set offset/scale values at the PAM level.
676 2 : poBand->SetOffset(CPLAtofM(poDS->GetKeyword("QUBE.CORE_BASE", "0.0")));
677 4 : poBand->SetScale(
678 2 : CPLAtofM(poDS->GetKeyword("QUBE.CORE_MULTIPLIER", "1.0")));
679 :
680 2 : poDS->SetBand(i + 1, std::move(poBand));
681 : }
682 :
683 : /* -------------------------------------------------------------------- */
684 : /* Check for a .prj file. For isis2 I would like to keep this in */
685 : /* -------------------------------------------------------------------- */
686 4 : const CPLString osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
687 4 : const CPLString osName = CPLGetBasenameSafe(poOpenInfo->pszFilename);
688 4 : const std::string osPrjFile = CPLFormCIFilenameSafe(osPath, osName, "prj");
689 :
690 2 : VSILFILE *fp = VSIFOpenL(osPrjFile.c_str(), "r");
691 2 : if (fp != nullptr)
692 : {
693 0 : VSIFCloseL(fp);
694 :
695 0 : char **papszLines = CSLLoad(osPrjFile.c_str());
696 :
697 0 : poDS->m_oSRS.importFromESRI(papszLines);
698 :
699 0 : CSLDestroy(papszLines);
700 : }
701 :
702 2 : if (dfULXMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0)
703 : {
704 2 : poDS->bGotTransform = TRUE;
705 2 : poDS->adfGeoTransform[0] = dfULXMap;
706 2 : poDS->adfGeoTransform[1] = dfXDim;
707 2 : poDS->adfGeoTransform[2] = 0.0;
708 2 : poDS->adfGeoTransform[3] = dfULYMap;
709 2 : poDS->adfGeoTransform[4] = 0.0;
710 2 : poDS->adfGeoTransform[5] = dfYDim;
711 : }
712 :
713 2 : if (!poDS->bGotTransform)
714 0 : poDS->bGotTransform = GDALReadWorldFile(poOpenInfo->pszFilename, "cbw",
715 0 : poDS->adfGeoTransform);
716 :
717 2 : if (!poDS->bGotTransform)
718 0 : poDS->bGotTransform = GDALReadWorldFile(poOpenInfo->pszFilename, "wld",
719 0 : poDS->adfGeoTransform);
720 :
721 : /* -------------------------------------------------------------------- */
722 : /* Initialize any PAM information. */
723 : /* -------------------------------------------------------------------- */
724 2 : poDS->SetDescription(poOpenInfo->pszFilename);
725 2 : poDS->TryLoadXML();
726 :
727 : /* -------------------------------------------------------------------- */
728 : /* Check for overviews. */
729 : /* -------------------------------------------------------------------- */
730 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
731 :
732 2 : return poDS.release();
733 : }
734 :
735 : /************************************************************************/
736 : /* GetKeyword() */
737 : /************************************************************************/
738 :
739 38 : const char *ISIS2Dataset::GetKeyword(const char *pszPath,
740 : const char *pszDefault)
741 :
742 : {
743 38 : return oKeywords.GetKeyword(pszPath, pszDefault);
744 : }
745 :
746 : /************************************************************************/
747 : /* GetKeywordSub() */
748 : /************************************************************************/
749 :
750 12 : const char *ISIS2Dataset::GetKeywordSub(const char *pszPath, int iSubscript,
751 : const char *pszDefault)
752 :
753 : {
754 12 : const char *pszResult = oKeywords.GetKeyword(pszPath, nullptr);
755 :
756 12 : if (pszResult == nullptr)
757 0 : return pszDefault;
758 :
759 12 : if (pszResult[0] != '(')
760 0 : return pszDefault;
761 :
762 : char **papszTokens =
763 12 : CSLTokenizeString2(pszResult, "(,)", CSLT_HONOURSTRINGS);
764 :
765 12 : if (iSubscript <= CSLCount(papszTokens))
766 : {
767 12 : oTempResult = papszTokens[iSubscript - 1];
768 12 : CSLDestroy(papszTokens);
769 12 : return oTempResult.c_str();
770 : }
771 :
772 0 : CSLDestroy(papszTokens);
773 0 : return pszDefault;
774 : }
775 :
776 : /************************************************************************/
777 : /* CleanString() */
778 : /* */
779 : /* Removes single or double quotes, and converts spaces to underscores. */
780 : /* The change is made in-place to CPLString. */
781 : /************************************************************************/
782 :
783 2 : void ISIS2Dataset::CleanString(CPLString &osInput)
784 :
785 : {
786 4 : if ((osInput.size() < 2) ||
787 2 : ((osInput.at(0) != '"' || osInput.back() != '"') &&
788 2 : (osInput.at(0) != '\'' || osInput.back() != '\'')))
789 2 : return;
790 :
791 0 : char *pszWrk = CPLStrdup(osInput.c_str() + 1);
792 :
793 0 : pszWrk[strlen(pszWrk) - 1] = '\0';
794 :
795 0 : for (int i = 0; pszWrk[i] != '\0'; i++)
796 : {
797 0 : if (pszWrk[i] == ' ')
798 0 : pszWrk[i] = '_';
799 : }
800 :
801 0 : osInput = pszWrk;
802 0 : CPLFree(pszWrk);
803 : }
804 :
805 : /************************************************************************/
806 : /* GDALRegister_ISIS2() */
807 : /************************************************************************/
808 :
809 1667 : void GDALRegister_ISIS2()
810 :
811 : {
812 1667 : if (GDALGetDriverByName(ISIS2_DRIVER_NAME) != nullptr)
813 282 : return;
814 :
815 1385 : GDALDriver *poDriver = new GDALDriver();
816 1385 : ISIS2DriverSetCommonMetadata(poDriver);
817 :
818 1385 : poDriver->pfnOpen = ISIS2Dataset::Open;
819 :
820 1385 : GetGDALDriverManager()->RegisterDriver(poDriver);
821 : }
|
