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