Merge pull request #1903 from rouault/vertical_perspective_projection

 ISIS3: add support for PointPerspective projection (refs #1856)

autotest/gcore/tiff_srs.py

@@ -592,3 +592,31 @@ def test_tiff_srs_read_epsg32631_4979_geotiff1_1():
     sr = ds.GetSpatialRef()
     # PROJ 6.0 didn't include the ID of the base CRS
     assert sr.ExportToWkt().replace(',ID["EPSG",4979]','') == 'PROJCRS["WGS 84 / UTM zone 31N",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4979]],CONVERSION["UTM zone 31N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",3,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]],ID["EPSG",16031]],CS[Cartesian,3],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1]]]'.replace(',ID["EPSG",4979]','')
+
+
+def test_tiff_srs_write_vertical_perspective():
+
+    ds = gdal.GetDriverByName('GTiff').Create('/vsimem/src.tif', 1, 1)
+    sr = osr.SpatialReference()
+    sr.SetGeogCS("GEOG_NAME", "D_DATUM_NAME", "", 3000000, 0)
+    sr.SetVerticalPerspective(1, 2, 0, 1000, 0, 0)
+    gdal.ErrorReset()
+    ds.SetSpatialRef(sr)
+    assert gdal.GetLastErrorMsg() == ''
+    ds = None
+
+    src_ds = gdal.Open('/vsimem/src.tif')
+    # First is PROJ 7
+    assert src_ds.GetSpatialRef().ExportToProj4() in ('+proj=nsper +lat_0=1 +lon_0=2 +h=1000 +x_0=0 +y_0=0 +R=3000000 +units=m +no_defs', '+proj=nsper +R=3000000 +lat_0=1 +lon_0=2 +h=1000 +x_0=0 +y_0=0 +wktext +no_defs')
+    gdal.ErrorReset()
+    gdal.GetDriverByName('GTiff').CreateCopy('/vsimem/dst.tif', src_ds)
+    assert gdal.GetLastErrorMsg() == ''
+
+    ds = gdal.Open('/vsimem/dst.tif')
+    assert ds.GetSpatialRef().ExportToProj4() == src_ds.GetSpatialRef().ExportToProj4()
+
+    src_ds = None
+    ds = None
+
+    gdal.GetDriverByName('GTiff').Delete('/vsimem/src.tif')
+    gdal.GetDriverByName('GTiff').Delete('/vsimem/dst.tif')

autotest/gdrivers/data/isis3_pointperspective.cub

@@ -0,0 +1,43 @@
+Object = IsisCube
+  Object = Core
+    Format      = Tile
+    TileSamples = 1244
+    TileLines   = 1244
+
+    Group = Dimensions
+      Samples = 1244
+      Lines   = 1244
+      Bands   = 3
+    End_Group
+
+    Group = Pixels
+      Type       = UnsignedByte
+      ByteOrder  = Lsb
+      Base       = 0.0
+      Multiplier = 1.0
+    End_Group
+  End_Object
+
+  Group = Mapping
+    ProjectionName     = PointPerspective
+    CenterLongitude    = -90.0
+    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   = -3110000.0 <meters>
+    UpperLeftCornerY   = 3110000.0 <meters>
+    PixelResolution    = 5000.0 <meters/pixel>
+    Scale              = 11.852768147075 <pixels/degree>
+    CenterLatitude     = -10.0
+    Distance           = 35000.0
+  End_Group
+End_Object
+
+End

autotest/gdrivers/isis.py

@@ -1687,3 +1687,32 @@ def test_isis3_preserve_label_across_format():
 
     src_ds = None
     gdal.Unlink('/vsimem/multiband.lbl')
+
+
+def test_isis3_point_perspective_read():
+
+    ds = gdal.Open('data/isis3_pointperspective.cub')
+    assert ds.GetSpatialRef().ExportToProj4() == '+proj=nsper +lat_0=-10 +lon_0=-90 +h=31603810 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs'
+
+
+def test_isis3_point_perspective_write():
+
+    if osr.GetPROJVersionMajor() < 7:
+        pytest.skip()
+
+    sr = osr.SpatialReference()
+    sr.SetGeogCS("GEOG_NAME", "D_DATUM_NAME", "", 3000000, 0)
+    sr.SetVerticalPerspective(1, 2, 0, 1000, 0, 0)
+    ds = gdal.GetDriverByName('ISIS3').Create('/vsimem/isis_tmp.lbl', 1, 1)
+    ds.SetSpatialRef(sr)
+    ds.SetGeoTransform([-10, 1, 0, 40, 0, -1])
+    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')

autotest/osr/osr_basic.py

@@ -1613,3 +1613,14 @@ def test_osr_promote_to_3D():
 
     assert sr.PromoteTo3D() == 0
     assert sr.GetAuthorityCode(None) == '4979'
+
+
+def test_osr_SetVerticalPerspective():
+
+    sr = osr.SpatialReference()
+    sr.SetVerticalPerspective(1, 2, 0, 3, 4, 5)
+    assert sr.ExportToProj4() == '+proj=nsper +lat_0=1 +lon_0=2 +h=3 +x_0=4 +y_0=5 +datum=WGS84 +units=m +no_defs'
+    if osr.GetPROJVersionMajor() >= 7:
+        assert sr.GetAttrValue('PROJECTION') in 'Vertical Perspective'
+        assert sr.GetNormProjParm('Longitude of topocentric origin') == 2
+

gdal/apps/gdal_translate_lib.cpp

@@ -691,7 +691,17 @@ GDALDatasetH GDALTranslate( const char *pszDest, GDALDatasetH hSrcDataset,
         }
 
         char* pszSRS = nullptr;
-        oOutputSRS.exportToWkt( &pszSRS );
+        {
+            CPLErrorStateBackuper oErrorStateBackuper;
+            CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
+            if( oOutputSRS.exportToWkt( &pszSRS ) != OGRERR_NONE )
+            {
+                CPLFree(pszSRS);
+                pszSRS = nullptr;
+                const char* const apszOptions[] = { "FORMAT=WKT2", nullptr };
+                oOutputSRS.exportToWkt( &pszSRS, apszOptions );
+            }
+        }
         CPLFree( psOptions->pszOutputSRS );
         psOptions->pszOutputSRS = CPLStrdup( pszSRS );
         CPLFree( pszSRS );

gdal/apps/gdalwarp_lib.cpp

@@ -290,22 +290,36 @@ static double GetAverageSegmentLength(OGRGeometryH hGeom)
 /* option to determine the source SRS.                                  */
 /************************************************************************/
 
-static const char* GetSrcDSProjection( GDALDatasetH hDS,
+static CPLString GetSrcDSProjection( GDALDatasetH hDS,
                                        char** papszTO )
 {
     const char *pszProjection = CSLFetchNameValue( papszTO, "SRC_SRS" );
     if( pszProjection != nullptr || hDS == nullptr )
     {
-        return pszProjection;
+        return pszProjection ? pszProjection : "";
     }
 
     const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" );
     char** papszMD = nullptr;
-    if( GDALGetProjectionRef( hDS ) != nullptr
-        && strlen(GDALGetProjectionRef( hDS )) > 0
+    const OGRSpatialReferenceH hSRS = GDALGetSpatialRef( hDS );
+    if( hSRS
         && (pszMethod == nullptr || EQUAL(pszMethod,"GEOTRANSFORM")) )
     {
-        pszProjection = GDALGetProjectionRef( hDS );
+        char* pszWKT = nullptr;
+        {
+            CPLErrorStateBackuper oErrorStateBackuper;
+            CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
+            if( OSRExportToWkt(hSRS, &pszWKT) != OGRERR_NONE )
+            {
+                CPLFree(pszWKT);
+                pszWKT = nullptr;
+                const char* const apszOptions[] = { "FORMAT=WKT2", nullptr };
+                OSRExportToWktEx(hSRS, &pszWKT, apszOptions );
+            }
+        }
+        CPLString osWKT = pszWKT ? pszWKT : "";
+        CPLFree(pszWKT);
+        return osWKT;
     }
     else if( GDALGetGCPProjection( hDS ) != nullptr
              && strlen(GDALGetGCPProjection( hDS )) > 0
@@ -324,7 +338,7 @@ static const char* GetSrcDSProjection( GDALDatasetH hDS,
     {
         pszProjection = CSLFetchNameValue( papszMD, "SRS" );
     }
-    return pszProjection;
+    return pszProjection ? pszProjection : "";
 }
 
 /************************************************************************/
@@ -348,15 +362,15 @@ static CPLErr CropToCutline( OGRGeometryH hCutline, char** papszTO,
     OGRSpatialReferenceH hSrcSRS = nullptr;
     OGRSpatialReferenceH hDstSRS = nullptr;
 
-    const char *pszThisSourceSRS =
+    const CPLString osThisSourceSRS =
         GetSrcDSProjection(
             nSrcCount > 0 ? pahSrcDS[0] : nullptr,
             papszTO);
-    if( pszThisSourceSRS != nullptr && pszThisSourceSRS[0] != '\0' )
+    if( !osThisSourceSRS.empty() )
     {
         hSrcSRS = OSRNewSpatialReference(nullptr);
         OSRSetAxisMappingStrategy(hSrcSRS, OAMS_TRADITIONAL_GIS_ORDER);
-        if( OSRSetFromUserInput( hSrcSRS, pszThisSourceSRS ) != OGRERR_NONE )
+        if( OSRSetFromUserInput( hSrcSRS, osThisSourceSRS ) != OGRERR_NONE )
         {
             CPLError(CE_Failure, CPLE_AppDefined,
                      "Cannot compute bounding box of cutline.");
@@ -553,22 +567,29 @@ static GDALDatasetH ApplyVerticalShiftGrid( GDALDatasetH hWrkSrcDS,
 {
     bErrorOccurredOut = false;
     // Check if we must do vertical shift grid transform
-    const char* pszSrcWKT = CSLFetchNameValueDef(
-                                    psOptions->papszTO, "SRC_SRS",
-                                    GDALGetProjectionRef(hWrkSrcDS) );
+    OGRSpatialReference oSRSSrc;
+    OGRSpatialReference oSRSDst;
+    const char* pszSrcWKT = CSLFetchNameValue(psOptions->papszTO, "SRC_SRS");
+    if( pszSrcWKT )
+        oSRSSrc.SetFromUserInput( pszSrcWKT );
+    else
+    {
+        auto hSRS = GDALGetSpatialRef(hWrkSrcDS);
+        if( hSRS )
+            oSRSSrc = *(OGRSpatialReference::FromHandle(hSRS));
+    }
+
     const char* pszDstWKT = CSLFetchNameValue( psOptions->papszTO, "DST_SRS" );
+    if( pszDstWKT )
+        oSRSDst.SetFromUserInput( pszDstWKT );
+
     double adfGT[6] = {};
     if( GDALGetRasterCount(hWrkSrcDS) == 1 &&
         GDALGetGeoTransform(hWrkSrcDS, adfGT) == CE_None &&
-        pszSrcWKT != nullptr && pszSrcWKT[0] != '\0' &&
-        pszDstWKT != nullptr )
+        !oSRSSrc.IsEmpty() && !oSRSDst.IsEmpty() )
     {
-        OGRSpatialReference oSRSSrc;
-        OGRSpatialReference oSRSDst;
-        if( oSRSSrc.SetFromUserInput( pszSrcWKT ) == OGRERR_NONE &&
-            oSRSDst.SetFromUserInput( pszDstWKT ) == OGRERR_NONE &&
-            (oSRSSrc.IsCompound() || Is3DGeogcs(oSRSSrc) ||
-             oSRSDst.IsCompound() || Is3DGeogcs(oSRSDst)) )
+        if( (oSRSSrc.IsCompound() || Is3DGeogcs(oSRSSrc)) ||
+            (oSRSDst.IsCompound() || Is3DGeogcs(oSRSDst)) )
         {
             const char *pszSrcProj4Geoids =
                 oSRSSrc.GetExtension( "VERT_DATUM", "PROJ4_GRIDS" );
@@ -3004,12 +3025,11 @@ GDALWarpCreateOutput( int nSrcCount, GDALDatasetH *pahSrcDS, const char *pszFile
 /* -------------------------------------------------------------------- */
         if( iSrc == 0 && osThisTargetSRS.empty() )
         {
-            const char* pszThisSourceSRS = GetSrcDSProjection(
-                                                    pahSrcDS[iSrc], papszTO );
-            if( pszThisSourceSRS != nullptr && pszThisSourceSRS[0] != '\0' )
+            const auto osThisSourceSRS = GetSrcDSProjection( pahSrcDS[iSrc], papszTO );
+            if( !osThisSourceSRS.empty() )
             {
-                osThisTargetSRS = pszThisSourceSRS;
-                aoTOList.SetNameValue("DST_SRS", pszThisSourceSRS);
+                osThisTargetSRS = osThisSourceSRS;
+                aoTOList.SetNameValue("DST_SRS", osThisSourceSRS);
             }
         }
 
@@ -3518,12 +3538,12 @@ TransformCutlineToSource( GDALDatasetH hSrcDS, OGRGeometryH hCutline,
 /*      one.                                                            */
 /* -------------------------------------------------------------------- */
     OGRSpatialReferenceH  hRasterSRS = nullptr;
-    const char *pszProjection = GetSrcDSProjection( hSrcDS, papszTO_In);
-    if( pszProjection != nullptr )
+    const CPLString osProjection = GetSrcDSProjection( hSrcDS, papszTO_In);
+    if( !osProjection.empty() )
     {
         hRasterSRS = OSRNewSpatialReference(nullptr);
         OSRSetAxisMappingStrategy(hRasterSRS, OAMS_TRADITIONAL_GIS_ORDER);
-        if( OSRSetFromUserInput( hRasterSRS, pszProjection ) != OGRERR_NONE )
+        if( OSRSetFromUserInput( hRasterSRS, osProjection ) != OGRERR_NONE )
         {
             OSRDestroySpatialReference(hRasterSRS);
             hRasterSRS = nullptr;

gdal/frmts/gtiff/geotiff.cpp

@@ -11123,11 +11123,21 @@ void GTiffDataset::WriteGeoTIFFInfo()
         if( bHasProjection )
         {
             char* pszProjection = nullptr;
-            m_oSRS.exportToWkt(&pszProjection);
+            {
+                CPLErrorStateBackuper oErrorStateBackuper;
+                CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
+                m_oSRS.exportToWkt(&pszProjection);
+            }
             if( pszProjection && pszProjection[0] )
+            {
                 GTIFSetFromOGISDefnEx( psGTIF, pszProjection,
                                        m_eGeoTIFFKeysFlavor,
                                        m_eGeoTIFFVersion );
+            }
+            else
+            {
+                GDALPamDataset::SetSpatialRef(&m_oSRS);
+            }
             CPLFree(pszProjection);
         }
 
@@ -17560,17 +17570,30 @@ GTiffDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
 /*      Write the projection information, if possible.                  */
 /* -------------------------------------------------------------------- */
     const bool bHasProjection = l_poSRS != nullptr;
+    bool bExportSRSToPAM = false;
     if( (bHasProjection || bPixelIsPoint) && bGeoTIFF )
     {
         GTIF *psGTIF = GTiffDatasetGTIFNew( l_hTIFF );
 
         if( bHasProjection )
         {
             char* pszWKT = nullptr;
-            l_poSRS->exportToWkt(&pszWKT);
-            GTIFSetFromOGISDefnEx( psGTIF, pszWKT,
-                                   GetGTIFFKeysFlavor(papszOptions),
-                                   GetGeoTIFFVersion(papszOptions) );
+            OGRErr eErr;
+            {
+                CPLErrorStateBackuper oErrorStateBackuper;
+                CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
+                eErr = l_poSRS->exportToWkt(&pszWKT);
+            }
+            if( eErr == OGRERR_NONE )
+            {
+                GTIFSetFromOGISDefnEx( psGTIF, pszWKT,
+                                    GetGTIFFKeysFlavor(papszOptions),
+                                    GetGeoTIFFVersion(papszOptions) );
+            }
+            else
+            {
+                bExportSRSToPAM = true;
+            }
             CPLFree(pszWKT);
         }
 
@@ -17759,7 +17782,7 @@ GTiffDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
         osOldGTIFF_REPORT_COMPD_CSVal.empty() ? nullptr :
         osOldGTIFF_REPORT_COMPD_CSVal.c_str());
 
-    if( !bGeoTIFF && (poDS->GetPamFlags() & GPF_DISABLED) == 0 )
+    if( (!bGeoTIFF || bExportSRSToPAM) && (poDS->GetPamFlags() & GPF_DISABLED) == 0 )
     {
         // Copy georeferencing info to PAM if the profile is not GeoTIFF
         poDS->GDALPamDataset::SetSpatialRef(poDS->GetSpatialRef());