ISIS3: add support for Oblique Cylindrical projection (refs OSGeo#1856)

autotest/gdrivers/data/isis3_obliquecylindrical.cub

@@ -0,0 +1,43 @@
+Object = IsisCube
+  Object = Core
+    Format      = Tile
+    TileSamples = 1067
+    TileLines   = 1067
+
+    Group = Dimensions
+      Samples = 3201
+      Lines   = 2134
+      Bands   = 1
+    End_Group
+
+    Group = Pixels
+      Type       = UnsignedByte
+      ByteOrder  = Lsb
+      Base       = 0.0
+      Multiplier = 1.0
+    End_Group
+  End_Object
+
+  Group = Mapping
+    ProjectionName     = ObliqueCylindrical
+    TargetName         = Mars
+    EquatorialRadius   = 3396190.0 <meters>
+    PolarRadius        = 3376200.0 <meters>
+    LatitudeType       = Planetocentric
+    LongitudeDirection = PositiveEast
+    LongitudeDomain    = 180
+    MinimumLatitude    = -90.0
+    MaximumLatitude    = 90.0
+    MinimumLongitude   = -180.0
+    MaximumLongitude   = 180.0
+    UpperLeftCornerX   = -10670000.0 <meters>
+    UpperLeftCornerY   = 5335000.0 <meters>
+    PixelResolution    = 5000.0 <meters/pixel>
+    Scale              = 11.854939504661 <pixels/degree>
+    PoleLatitude       = 0.0
+    PoleLongitude      = 0.0
+    PoleRotation       = 90.0
+  End_Group
+End_Object
+
+End

autotest/gdrivers/isis.py

@@ -1709,10 +1709,46 @@ def test_isis3_point_perspective_write():
     ds = None
     ds = gdal.Open('/vsimem/isis_tmp.lbl')
     lbl = ds.GetMetadata_List('json:ISIS3')[0]
-    print(lbl)
     assert '"CenterLatitude":1.0' in lbl
     assert '"CenterLongitude":2.0' in lbl
     assert '"Distance":3001.0' in lbl
     ds = None
 
     gdal.GetDriverByName('ISIS3').Delete('/vsimem/isis_tmp.lbl')
+
+
+def test_isis3_oblique_cylindrical_read():
+
+    ds = gdal.Open('data/isis3_obliquecylindrical.cub')
+    srs = ds.GetSpatialRef()
+    assert srs.ExportToProj4() == '+proj=ob_tran +o_proj=eqc +o_lon_p=-90 +o_lat_p=180 +lon_0=0 +R=3396190 +units=m +no_defs'
+
+    pixel = ds.RasterXSize / 2.0
+    line = ds.RasterYSize / 2.0
+    gt = ds.GetGeoTransform()
+    x = gt[0] + pixel * gt[1] + line * gt[2]
+    y = gt[3] + pixel * gt[4] + line * gt[5]
+    geog_srs = srs.CloneGeogCS()
+    ct = osr.CoordinateTransformation(srs, geog_srs)
+    lon, lat, _ = ct.TransformPoint(x, y)
+    assert lon == pytest.approx(90.0)
+    assert lat == pytest.approx(-45.0, 1e-2)
+
+
+def test_isis3_oblique_cylindrical_write():
+
+    src_ds = gdal.Open('data/isis3_obliquecylindrical.cub')
+    ds = gdal.GetDriverByName('ISIS3').Create('/vsimem/isis_tmp.lbl', src_ds.RasterXSize, src_ds.RasterYSize)
+    ds.SetSpatialRef(src_ds.GetSpatialRef())
+    ds.SetGeoTransform(src_ds.GetGeoTransform())
+    src_ds = None
+    ds = None
+
+    ds = gdal.Open('/vsimem/isis_tmp.lbl')
+    lbl = ds.GetMetadata_List('json:ISIS3')[0]
+    assert '"PoleLongitude":0.0' in lbl
+    assert '"PoleLatitude":0.0' in lbl
+    assert '"PoleRotation":90.0' in lbl
+    ds = None
+
+    gdal.GetDriverByName('ISIS3').Delete('/vsimem/isis_tmp.lbl')

gdal/frmts/pds/isis3dataset.cpp

@@ -1880,6 +1880,13 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo )
     /**** This is the planets name i.e. Mars ***/
     const char *target_name = poDS->GetKeyword("IsisCube.Mapping.TargetName");
 
+#ifdef notdef
+    const double dfLongitudeMulFactor =
+        EQUAL(poDS->GetKeyword( "IsisCube.Mapping.LongitudeDirection", "PositiveEast"), "PositiveEast") ? 1 : -1;
+#else
+    const double dfLongitudeMulFactor = 1;
+#endif
+
     /***********   Grab MAP_PROJECTION_TYPE ************/
      const char *map_proj_name =
         poDS->GetKeyword( "IsisCube.Mapping.ProjectionName");
@@ -1898,7 +1905,7 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo )
 
     /***********   Grab CENTER_LON ************/
     const double center_lon =
-        CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.CenterLongitude"));
+        CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.CenterLongitude")) * dfLongitudeMulFactor;
 
     /***********   Grab 1st std parallel ************/
     const double first_std_parallel =
@@ -1929,7 +1936,7 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo )
     //  Equirectangular
     //  LambertConformal
     //  Mercator
-    //  ObliqueCylindrical //Todo
+    //  ObliqueCylindrical
     //  Orthographic
     //  PolarStereographic
     //  SimpleCylindrical
@@ -1959,11 +1966,28 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo )
     } else if (EQUAL( map_proj_name, "LambertConformal" )) {
         oSRS.SetLCC ( first_std_parallel, second_std_parallel, center_lat, center_lon, 0, 0 );
     } else if (EQUAL( map_proj_name, "PointPerspective" )) {
-        CPLString oProj4String;
         // Distance parameter is the distance to the center of the body, and is given in km
         const double distance = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.Distance")) * 1000.0;
         const double height_above_ground = distance - semi_major;
         oSRS.SetVerticalPerspective(center_lat, center_lon, 0, height_above_ground, 0, 0);
+    } else if (EQUAL( map_proj_name, "ObliqueCylindrical" )) {
+        const double poleLatitude = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PoleLatitude"));
+        const double poleLongitude = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PoleLongitude")) * dfLongitudeMulFactor;
+        const double poleRotation = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PoleRotation"));
+        CPLString oProj4String;
+        // ISIS3 rotated pole doesn't use the same conventions than PROJ ob_tran
+        // Compare the sign difference in https://github.com/USGS-Astrogeology/ISIS3/blob/3.8.0/isis/src/base/objs/ObliqueCylindrical/ObliqueCylindrical.cpp#L244
+        // and https://github.com/OSGeo/PROJ/blob/6.2/src/projections/ob_tran.cpp#L34
+        // They can be compensated by modifying the poleLatitude to 180-poleLatitude
+        // There's also a sign difference for the poleRotation parameter
+        // The existence of those different conventions is acknowledged in
+        // https://pds-imaging.jpl.nasa.gov/documentation/Cassini_BIDRSIS.PDF in the middle of page 10
+        oProj4String.Printf(
+            "+proj=ob_tran +o_proj=eqc +o_lon_p=%.18g +o_lat_p=%.18g +lon_0=%.18g",
+            -poleRotation,
+            180-poleLatitude,
+            poleLongitude);
+        oSRS.SetFromUserInput(oProj4String);
     } else {
         CPLDebug( "ISIS3",
                   "Dataset projection %s is not supported. Continuing...",
@@ -2914,6 +2938,40 @@ void ISIS3Dataset::BuildLabel()
                         (oSRS.GetNormProjParm("Viewpoint height", 0.0) + oSRS.GetSemiMajor()) / 1000.0 );
             }
 
+            else if( EQUAL(pszProjection, "custom_proj4") )
+            {
+                const char* pszProj4 = oSRS.GetExtension("PROJCS", "PROJ4", nullptr);
+                if( pszProj4 && strstr(pszProj4, "+proj=ob_tran" ) &&
+                    strstr(pszProj4, "+o_proj=eqc") )
+                {
+                    const auto FetchParam = [](const char* pszProj4Str, const char* pszKey)
+                    {
+                        CPLString needle;
+                        needle.Printf("+%s=", pszKey);
+                        const char* pszVal = strstr(pszProj4Str, needle.c_str());
+                        if( pszVal )
+                            return CPLAtof(pszVal+needle.size());
+                        return 0.0;
+                    };
+
+                    double dfLonP = FetchParam(pszProj4, "o_lon_p");
+                    double dfLatP = FetchParam(pszProj4, "o_lat_p");
+                    double dfLon0 = FetchParam(pszProj4, "lon_0");
+                    double dfPoleRotation = -dfLonP;
+                    double dfPoleLatitude = 180 - dfLatP;
+                    double dfPoleLongitude = dfLon0;
+                    oMapping.Add( "ProjectionName", "ObliqueCylindrical" );
+                    oMapping.Add( "PoleLatitude", dfPoleLatitude );
+                    oMapping.Add( "PoleLongitude", FixLong(dfPoleLongitude) );
+                    oMapping.Add( "PoleRotation", dfPoleRotation );
+                }
+                else
+                {
+                    CPLError(CE_Warning, CPLE_NotSupported,
+                            "Projection %s not supported",
+                            pszProjection);
+                }
+            }
             else
             {
                 CPLError(CE_Warning, CPLE_NotSupported,