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