Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: OGRSpatialReference translation from an ISIS3 (PVL) Mapping group.
5 : * Factored out of the ISIS3 driver (frmts/pds/isis3dataset.cpp).
6 : * Author: Trent Hare (thare@usgs.gov)
7 : * Frank Warmerdam (warmerdam@pobox.com)
8 : * Even Rouault (even.rouault at spatialys.com)
9 : * Oleg Alexandrov
10 : *
11 : * NOTE: Original code authored by Trent and placed in the public domain as
12 : * per US government policy. I have (within my rights) appropriated it and
13 : * placed it under the following license. This is not intended to diminish
14 : * Trents contribution.
15 : ******************************************************************************
16 : * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
17 : * Copyright (c) 2009-2026, Even Rouault <even.rouault at spatialys.com>
18 : * Copyright (c) 2017 Hobu Inc
19 : * Copyright (c) 2017, Dmitry Baryshnikov <polimax@mail.ru>
20 : * Copyright (c) 2017, NextGIS <info@nextgis.com>
21 : * Copyright (c) 2026, Oleg Alexandrov
22 : *
23 : * SPDX-License-Identifier: MIT
24 : ****************************************************************************/
25 :
26 : #include "cpl_port.h"
27 : #include "ogr_srs_api.h"
28 :
29 : #include <cmath>
30 : #include <string>
31 :
32 : #include "cpl_conv.h"
33 : #include "cpl_error.h"
34 : #include "cpl_json.h"
35 : #include "cpl_string.h"
36 : #include "ogr_spatialref.h"
37 :
38 : /************************************************************************/
39 : /* GetISISMappingString() */
40 : /************************************************************************/
41 :
42 : // Read a value from the ISIS PVL Mapping group as a string. A value that
43 : // carries a unit is stored as a { "value", "unit" } object, so unwrap that.
44 :
45 532 : static CPLString GetISISMappingString(const CPLJSONObject &oMapping,
46 : const char *pszKey,
47 : const char *pszDefault = "")
48 : {
49 1596 : CPLJSONObject oChild = oMapping.GetObj(pszKey);
50 532 : if (oChild.GetType() == CPLJSONObject::Type::Object)
51 0 : oChild = oChild.GetObj("value");
52 532 : if (oChild.GetType() == CPLJSONObject::Type::String)
53 776 : return oChild.ToString(pszDefault);
54 144 : return pszDefault;
55 : }
56 :
57 : /************************************************************************/
58 : /* GetISISMappingDouble() */
59 : /************************************************************************/
60 :
61 : // Read a value from the ISIS PVL Mapping group as a double.
62 :
63 948 : static double GetISISMappingDouble(const CPLJSONObject &oMapping,
64 : const char *pszKey, double dfDefault = 0.0)
65 : {
66 2844 : CPLJSONObject oChild = oMapping.GetObj(pszKey);
67 948 : if (oChild.GetType() == CPLJSONObject::Type::Object)
68 238 : oChild = oChild.GetObj("value");
69 948 : switch (oChild.GetType())
70 : {
71 0 : case CPLJSONObject::Type::Integer:
72 : case CPLJSONObject::Type::Long:
73 0 : return static_cast<double>(oChild.ToLong());
74 508 : case CPLJSONObject::Type::Double:
75 508 : return oChild.ToDouble();
76 46 : case CPLJSONObject::Type::String:
77 46 : return CPLAtof(oChild.ToString().c_str());
78 394 : default:
79 394 : return dfDefault;
80 : }
81 : }
82 :
83 : /************************************************************************/
84 : /* ParseISISPVLMappingGroup() */
85 : /************************************************************************/
86 :
87 : // Parse an ISIS PVL Mapping group as a flat CPLJSONObject of keys and string
88 : // values. The "Group = Mapping" and "End_Group" lines are ignored, and a
89 : // quoted value (such as ProjStr) may span lines.
90 :
91 12 : static CPLJSONObject ParseISISPVLMappingGroup(const char *pszText)
92 : {
93 12 : CPLJSONObject oMapping;
94 12 : const char *pszIter = pszText;
95 107 : while (*pszIter)
96 : {
97 288 : while (*pszIter && isspace(static_cast<unsigned char>(*pszIter)))
98 193 : ++pszIter;
99 95 : if (!*pszIter)
100 0 : break;
101 : // A comment line.
102 95 : if (*pszIter == '#')
103 : {
104 0 : while (*pszIter && *pszIter != '\n')
105 0 : ++pszIter;
106 24 : continue;
107 : }
108 :
109 : // Read the keyword name.
110 95 : std::string osKey;
111 1193 : while (*pszIter && !isspace(static_cast<unsigned char>(*pszIter)) &&
112 1098 : *pszIter != '=')
113 : {
114 1098 : osKey += *pszIter;
115 1098 : ++pszIter;
116 : }
117 326 : while (*pszIter && isspace(static_cast<unsigned char>(*pszIter)))
118 231 : ++pszIter;
119 :
120 : // Group/Object terminators carry no value.
121 178 : if (EQUAL(osKey.c_str(), "End") || EQUAL(osKey.c_str(), "End_Group") ||
122 83 : EQUAL(osKey.c_str(), "End_Object"))
123 12 : continue;
124 :
125 83 : if (*pszIter != '=')
126 : {
127 : // A keyword with no value; skip to the end of the line.
128 0 : while (*pszIter && *pszIter != '\n')
129 0 : ++pszIter;
130 0 : continue;
131 : }
132 83 : ++pszIter; // skip '='
133 166 : while (*pszIter && (*pszIter == ' ' || *pszIter == '\t'))
134 83 : ++pszIter;
135 :
136 : // Read the value.
137 83 : std::string osValue;
138 83 : if (*pszIter == '"')
139 : {
140 2 : ++pszIter; // opening quote
141 158 : while (*pszIter && *pszIter != '"')
142 : {
143 156 : osValue += *pszIter;
144 156 : ++pszIter;
145 : }
146 2 : if (*pszIter == '"')
147 2 : ++pszIter; // closing quote
148 : }
149 : else
150 : {
151 842 : while (*pszIter && *pszIter != '\n' && *pszIter != '\r')
152 : {
153 761 : osValue += *pszIter;
154 761 : ++pszIter;
155 : }
156 162 : while (!osValue.empty() &&
157 81 : (osValue.back() == ' ' || osValue.back() == '\t'))
158 0 : osValue.pop_back();
159 : }
160 :
161 : // The "Group = Mapping" / "Object = ..." preamble is not a real key.
162 83 : if (EQUAL(osKey.c_str(), "Group") || EQUAL(osKey.c_str(), "Object"))
163 12 : continue;
164 :
165 71 : if (!osKey.empty())
166 71 : oMapping.Add(osKey, osValue);
167 : }
168 12 : return oMapping;
169 : }
170 :
171 : /************************************************************************/
172 : /* importFromISISPVL() */
173 : /************************************************************************/
174 :
175 : /**
176 : * \brief Import a coordinate system from an ISIS3 PVL Mapping group.
177 : *
178 : * This method reads the projection of an ISIS3 cube from the text of its
179 : * Mapping group. It supports the named ISIS projections as well as the newer
180 : * generic form where ProjectionName is IProj and the definition is carried by
181 : * a ProjStr PROJ string.
182 : *
183 : * This method is the equivalent of the C function OSRImportFromISISPVL().
184 : *
185 : * @param pszPVLMappingGroup text of the ISIS PVL Mapping group.
186 : *
187 : * @return OGRERR_NONE on success, or an error code on failure.
188 : *
189 : * @since GDAL 3.14
190 : */
191 :
192 12 : OGRErr OGRSpatialReference::importFromISISPVL(const char *pszPVLMappingGroup)
193 : {
194 12 : return importFromISISPVL(ParseISISPVLMappingGroup(pszPVLMappingGroup));
195 : }
196 :
197 : /************************************************************************/
198 : /* importFromISISPVL() */
199 : /************************************************************************/
200 :
201 : /**
202 : * \brief Import a coordinate system from an already parsed ISIS3 Mapping group.
203 : *
204 : * This is the overload of importFromISISPVL(const char*) that takes the
205 : * Mapping group as a CPLJSONObject of key/value pairs, as parsed by the
206 : * ISIS3 driver. See importFromISISPVL(const char*) for details.
207 : *
208 : * @param oMapping the ISIS Mapping group as a CPLJSONObject.
209 : *
210 : * @return OGRERR_NONE on success, or an error code on failure.
211 : *
212 : * @since GDAL 3.14
213 : */
214 :
215 133 : OGRErr OGRSpatialReference::importFromISISPVL(const CPLJSONObject &oMapping)
216 : {
217 133 : Clear();
218 :
219 : /*********** Grab TARGET_NAME ************/
220 : /**** This is the planets name i.e. Mars ***/
221 266 : const CPLString target_name = GetISISMappingString(oMapping, "TargetName");
222 :
223 : #ifdef notdef
224 : const double dfLongitudeMulFactor =
225 : EQUAL(GetISISMappingString(oMapping, "LongitudeDirection",
226 : "PositiveEast"),
227 : "PositiveEast")
228 : ? 1
229 : : -1;
230 : #else
231 133 : const double dfLongitudeMulFactor = 1;
232 : #endif
233 :
234 : /*********** Grab MAP_PROJECTION_TYPE ************/
235 : const CPLString map_proj_name =
236 266 : GetISISMappingString(oMapping, "ProjectionName");
237 :
238 : /*********** Grab SEMI-MAJOR ************/
239 : const double semi_major =
240 133 : GetISISMappingDouble(oMapping, "EquatorialRadius");
241 :
242 : /*********** Grab semi-minor ************/
243 133 : const double semi_minor = GetISISMappingDouble(oMapping, "PolarRadius");
244 :
245 : /*********** Grab CENTER_LAT ************/
246 133 : const double center_lat = GetISISMappingDouble(oMapping, "CenterLatitude");
247 :
248 : /*********** Grab CENTER_LON ************/
249 : const double center_lon =
250 133 : GetISISMappingDouble(oMapping, "CenterLongitude") *
251 : dfLongitudeMulFactor;
252 :
253 : /*********** Grab 1st std parallel ************/
254 : const double first_std_parallel =
255 133 : GetISISMappingDouble(oMapping, "FirstStandardParallel");
256 :
257 : /*********** Grab 2nd std parallel ************/
258 : const double second_std_parallel =
259 133 : GetISISMappingDouble(oMapping, "SecondStandardParallel");
260 :
261 : /*********** Grab scaleFactor ************/
262 : const double scaleFactor =
263 133 : GetISISMappingDouble(oMapping, "scaleFactor", 1.0);
264 :
265 : /*** grab LatitudeType = Planetographic ****/
266 : // Need to further study how ocentric/ographic will effect the gdal library
267 : // So far we will use this fact to define a sphere or ellipse for some
268 : // projections
269 :
270 : // Frank - may need to talk this over
271 133 : bool bIsGeographic = true;
272 133 : if (EQUAL(GetISISMappingString(oMapping, "LatitudeType"), "Planetocentric"))
273 92 : bIsGeographic = false;
274 :
275 : // Set oSRS projection and parameters
276 : // ############################################################
277 : // ISIS3 Projection types
278 : // Equirectangular
279 : // LambertConformal
280 : // Mercator
281 : // ObliqueCylindrical
282 : // Orthographic
283 : // PolarStereographic
284 : // SimpleCylindrical
285 : // Sinusoidal
286 : // TransverseMercator
287 :
288 : #ifdef DEBUG
289 133 : CPLDebug("ISIS3", "using projection %s", map_proj_name.c_str());
290 : #endif
291 :
292 133 : bool bProjectionSet = true;
293 :
294 : // The ISIS IProj projection (ProjectionName = IProj) carries its full
295 : // definition in a ProjStr PROJ string. When present, use it directly and
296 : // skip the per-parameter reconstruction below.
297 266 : const CPLString osProjStrRaw = GetISISMappingString(oMapping, "ProjStr");
298 133 : if (!osProjStrRaw.empty())
299 : {
300 : // A ProjStr that spans several label lines comes back with the line
301 : // breaks embedded as literal "\n" (plus indentation), or as real
302 : // control characters. Turn those into spaces before handing it to PROJ.
303 10 : std::string osProjStr;
304 444 : for (const char *pszIter = osProjStrRaw.c_str(); *pszIter; ++pszIter)
305 : {
306 439 : if (pszIter[0] == '\\' &&
307 2 : (pszIter[1] == 'n' || pszIter[1] == 'r' || pszIter[1] == 't'))
308 : {
309 2 : osProjStr += ' ';
310 2 : ++pszIter;
311 : }
312 437 : else if (*pszIter == '\n' || *pszIter == '\r' || *pszIter == '\t')
313 1 : osProjStr += ' ';
314 : else
315 436 : osProjStr += *pszIter;
316 : }
317 :
318 : OGRErr eErr;
319 : {
320 : // A ProjStr that PROJ rejects should leave a warning and an error
321 : // code, not a hard error that callers see as an exception.
322 5 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
323 5 : eErr = importFromProj4(osProjStr.c_str());
324 : }
325 5 : if (eErr == OGRERR_NONE)
326 : {
327 3 : SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
328 3 : return OGRERR_NONE;
329 : }
330 :
331 2 : CPLError(CE_Warning, CPLE_AppDefined,
332 : "Cannot parse ProjStr '%s' from ISIS3 Mapping group.",
333 : osProjStr.c_str());
334 2 : return OGRERR_CORRUPT_DATA;
335 : }
336 188 : else if ((EQUAL(map_proj_name, "Equirectangular")) ||
337 60 : (EQUAL(map_proj_name, "SimpleCylindrical")))
338 : {
339 84 : SetEquirectangular2(0.0, center_lon, center_lat, 0, 0);
340 : }
341 44 : else if (EQUAL(map_proj_name, "Orthographic"))
342 : {
343 3 : SetOrthographic(center_lat, center_lon, 0, 0);
344 : }
345 41 : else if (EQUAL(map_proj_name, "Sinusoidal"))
346 : {
347 3 : SetSinusoidal(center_lon, 0, 0);
348 : }
349 38 : else if (EQUAL(map_proj_name, "Mercator"))
350 : {
351 3 : SetMercator(center_lat, center_lon, scaleFactor, 0, 0);
352 : }
353 35 : else if (EQUAL(map_proj_name, "PolarStereographic"))
354 : {
355 3 : SetPS(center_lat, center_lon, scaleFactor, 0, 0);
356 : }
357 32 : else if (EQUAL(map_proj_name, "TransverseMercator"))
358 : {
359 22 : SetTM(center_lat, center_lon, scaleFactor, 0, 0);
360 : }
361 10 : else if (EQUAL(map_proj_name, "LambertConformal"))
362 : {
363 3 : SetLCC(first_std_parallel, second_std_parallel, center_lat, center_lon,
364 : 0, 0);
365 : }
366 7 : else if (EQUAL(map_proj_name, "PointPerspective"))
367 : {
368 : // Distance parameter is the distance to the center of the body, and is
369 : // given in km
370 : const double distance =
371 2 : GetISISMappingDouble(oMapping, "Distance") * 1000.0;
372 2 : const double height_above_ground = distance - semi_major;
373 2 : SetVerticalPerspective(center_lat, center_lon, 0, height_above_ground,
374 : 0, 0);
375 : }
376 5 : else if (EQUAL(map_proj_name, "ObliqueCylindrical"))
377 : {
378 : const double poleLatitude =
379 5 : GetISISMappingDouble(oMapping, "PoleLatitude");
380 : const double poleLongitude =
381 5 : GetISISMappingDouble(oMapping, "PoleLongitude") *
382 : dfLongitudeMulFactor;
383 : const double poleRotation =
384 5 : GetISISMappingDouble(oMapping, "PoleRotation");
385 10 : CPLString oProj4String;
386 : // ISIS3 rotated pole doesn't use the same conventions than PROJ ob_tran
387 : // Compare the sign difference in
388 : // https://github.com/USGS-Astrogeology/ISIS3/blob/3.8.0/isis/src/base/objs/ObliqueCylindrical/ObliqueCylindrical.cpp#L244
389 : // and
390 : // https://github.com/OSGeo/PROJ/blob/6.2/src/projections/ob_tran.cpp#L34
391 : // They can be compensated by modifying the poleLatitude to
392 : // 180-poleLatitude There's also a sign difference for the poleRotation
393 : // parameter The existence of those different conventions is
394 : // acknowledged in
395 : // https://pds-imaging.jpl.nasa.gov/documentation/Cassini_BIDRSIS.PDF in
396 : // the middle of page 10
397 : oProj4String.Printf("+proj=ob_tran +o_proj=eqc +o_lon_p=%.17g "
398 : "+o_lat_p=%.17g +lon_0=%.17g",
399 5 : -poleRotation, 180 - poleLatitude, poleLongitude);
400 5 : SetFromUserInput(oProj4String);
401 : }
402 : else
403 : {
404 0 : CPLDebug("ISIS3",
405 : "Dataset projection %s is not supported. Continuing...",
406 : map_proj_name.c_str());
407 0 : bProjectionSet = false;
408 : }
409 :
410 128 : if (!bProjectionSet)
411 0 : return OGRERR_UNSUPPORTED_SRS;
412 :
413 : // Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
414 256 : CPLString osProjTargetName(map_proj_name);
415 128 : osProjTargetName += " ";
416 128 : osProjTargetName += target_name;
417 128 : SetProjCS(osProjTargetName); // set ProjCS keyword
418 :
419 : // The geographic/geocentric name will be the same basic name as the
420 : // body name 'GCS' = Geographic/Geocentric Coordinate System
421 256 : CPLString osGeogName("GCS_");
422 128 : osGeogName += target_name;
423 :
424 : // The datum name will be the same basic name as the planet
425 256 : CPLString osDatumName("D_");
426 128 : osDatumName += target_name;
427 :
428 128 : CPLString osSphereName(target_name);
429 : // strcat(osSphereName, "_IAU_IAG"); //Might not be IAU defined so
430 : // don't add
431 :
432 : // calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
433 128 : double iflattening = 0.0;
434 128 : if ((semi_major - semi_minor) < 0.0000001)
435 38 : iflattening = 0;
436 : else
437 90 : iflattening = semi_major / (semi_major - semi_minor);
438 :
439 : // Set the body size but take into consideration which proj is being
440 : // used to help w/ proj4 compatibility The use of a Sphere, polar radius
441 : // or ellipse here is based on how ISIS does it internally
442 256 : if (((EQUAL(map_proj_name, "Stereographic") && (fabs(center_lat) == 90))) ||
443 128 : (EQUAL(map_proj_name, "PolarStereographic")))
444 : {
445 3 : if (bIsGeographic)
446 : {
447 : // Geograpraphic, so set an ellipse
448 1 : SetGeogCS(osGeogName, osDatumName, osSphereName, semi_major,
449 : iflattening, "Reference_Meridian", 0.0);
450 : }
451 : else
452 : {
453 : // Geocentric, so force a sphere using the semi-minor axis. I
454 : // hope...
455 2 : osSphereName += "_polarRadius";
456 2 : SetGeogCS(osGeogName, osDatumName, osSphereName, semi_minor, 0.0,
457 : "Reference_Meridian", 0.0);
458 : }
459 : }
460 125 : else if ((EQUAL(map_proj_name, "SimpleCylindrical")) ||
461 109 : (EQUAL(map_proj_name, "Orthographic")) ||
462 106 : (EQUAL(map_proj_name, "Stereographic")) ||
463 337 : (EQUAL(map_proj_name, "Sinusoidal")) ||
464 103 : (EQUAL(map_proj_name, "PointPerspective")))
465 : {
466 : // ISIS uses the spherical equation for these projections
467 : // so force a sphere.
468 24 : SetGeogCS(osGeogName, osDatumName, osSphereName, semi_major, 0.0,
469 : "Reference_Meridian", 0.0);
470 : }
471 101 : else if (EQUAL(map_proj_name, "Equirectangular"))
472 : {
473 : // Calculate localRadius using ISIS3 simple elliptical method
474 : // not the more standard Radius of Curvature method
475 : // PI = 4 * atan(1);
476 68 : const double radLat = center_lat * M_PI / 180; // in radians
477 68 : const double meanRadius = sqrt(pow(semi_minor * cos(radLat), 2) +
478 68 : pow(semi_major * sin(radLat), 2));
479 68 : const double localRadius =
480 68 : (meanRadius == 0.0) ? 0.0 : semi_major * semi_minor / meanRadius;
481 68 : osSphereName += "_localRadius";
482 68 : SetGeogCS(osGeogName, osDatumName, osSphereName, localRadius, 0.0,
483 : "Reference_Meridian", 0.0);
484 : }
485 : else
486 : {
487 : // All other projections: Mercator, Transverse Mercator, Lambert
488 : // Conformal, etc. Geographic, so set an ellipse
489 33 : if (bIsGeographic)
490 : {
491 4 : SetGeogCS(osGeogName, osDatumName, osSphereName, semi_major,
492 : iflattening, "Reference_Meridian", 0.0);
493 : }
494 : else
495 : {
496 : // Geocentric, so force a sphere. I hope...
497 29 : SetGeogCS(osGeogName, osDatumName, osSphereName, semi_major, 0.0,
498 : "Reference_Meridian", 0.0);
499 : }
500 : }
501 :
502 128 : SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
503 :
504 128 : return OGRERR_NONE;
505 : }
506 :
507 : /************************************************************************/
508 : /* OSRImportFromISISPVL() */
509 : /************************************************************************/
510 :
511 : /**
512 : * \brief Import a coordinate system from an ISIS3 PVL Mapping group.
513 : *
514 : * This function is the same as OGRSpatialReference::importFromISISPVL().
515 : *
516 : * @since GDAL 3.14
517 : */
518 :
519 12 : OGRErr OSRImportFromISISPVL(OGRSpatialReferenceH hSRS,
520 : const char *pszPVLMappingGroup)
521 : {
522 12 : VALIDATE_POINTER1(hSRS, "OSRImportFromISISPVL", OGRERR_FAILURE);
523 :
524 12 : return OGRSpatialReference::FromHandle(hSRS)->importFromISISPVL(
525 12 : pszPVLMappingGroup);
526 : }
|
