From e6f4c6697477ce23c3debabe21fae07d5c74c84c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 2 Oct 2019 22:59:28 +0200 Subject: [PATCH] ISIS3: add support for PointPerspective projection (refs #1856) --- .../gdrivers/data/isis3_pointperspective.cub | 43 ++++++++ autotest/gdrivers/isis.py | 29 ++++++ gdal/frmts/pds/isis3dataset.cpp | 97 +++++++++++-------- 3 files changed, 126 insertions(+), 43 deletions(-) create mode 100644 autotest/gdrivers/data/isis3_pointperspective.cub diff --git a/autotest/gdrivers/data/isis3_pointperspective.cub b/autotest/gdrivers/data/isis3_pointperspective.cub new file mode 100644 index 000000000000..bd4edf5add99 --- /dev/null +++ b/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 + PolarRadius = 3376200.0 + LatitudeType = Planetocentric + LongitudeDirection = PositiveEast + LongitudeDomain = 180 + MinimumLatitude = -90.0 + MaximumLatitude = 90.0 + MinimumLongitude = -180.0 + MaximumLongitude = 180.0 + UpperLeftCornerX = -3110000.0 + UpperLeftCornerY = 3110000.0 + PixelResolution = 5000.0 + Scale = 11.852768147075 + CenterLatitude = -10.0 + Distance = 35000.0 + End_Group +End_Object + +End \ No newline at end of file diff --git a/autotest/gdrivers/isis.py b/autotest/gdrivers/isis.py index 59b011df45b8..3f58a1d20fb9 100755 --- a/autotest/gdrivers/isis.py +++ b/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') diff --git a/gdal/frmts/pds/isis3dataset.cpp b/gdal/frmts/pds/isis3dataset.cpp index a08a13b1be06..6bf16ffdf0b7 100644 --- a/gdal/frmts/pds/isis3dataset.cpp +++ b/gdal/frmts/pds/isis3dataset.cpp @@ -151,7 +151,7 @@ class ISIS3Dataset final: public RawDataset bool m_bHasSrcNoData; // creation only double m_dfSrcNoData; // creation only - CPLString m_osProjection; + OGRSpatialReference m_oSRS; // creation only variables CPLString m_osComment; @@ -199,14 +199,8 @@ class ISIS3Dataset final: public RawDataset virtual CPLErr GetGeoTransform( double * padfTransform ) override; virtual CPLErr SetGeoTransform( double * padfTransform ) override; - virtual const char *_GetProjectionRef(void) override; - virtual CPLErr _SetProjection( const char* pszProjection ) override; - const OGRSpatialReference* GetSpatialRef() const override { - return GetSpatialRefFromOldGetProjectionRef(); - } - CPLErr SetSpatialRef(const OGRSpatialReference* poSRS) override { - return OldSetProjectionFromSetSpatialRef(poSRS); - } + const OGRSpatialReference* GetSpatialRef() const override; + CPLErr SetSpatialRef(const OGRSpatialReference* poSRS) override; virtual char **GetFileList() override; @@ -1441,29 +1435,32 @@ char **ISIS3Dataset::GetFileList() } /************************************************************************/ -/* GetProjectionRef() */ +/* GetSpatialRef() */ /************************************************************************/ -const char *ISIS3Dataset::_GetProjectionRef() +const OGRSpatialReference* ISIS3Dataset::GetSpatialRef() const { - if( !m_osProjection.empty() ) - return m_osProjection; + if( !m_oSRS.IsEmpty() ) + return &m_oSRS; - return GDALPamDataset::_GetProjectionRef(); + return GDALPamDataset::GetSpatialRef(); } /************************************************************************/ -/* SetProjection() */ +/* SetSpatialRef() */ /************************************************************************/ -CPLErr ISIS3Dataset::_SetProjection( const char* pszProjection ) +CPLErr ISIS3Dataset::SetSpatialRef( const OGRSpatialReference* poSRS ) { if( eAccess == GA_ReadOnly ) - return GDALPamDataset::_SetProjection( pszProjection ); - m_osProjection = pszProjection ? pszProjection : ""; + return GDALPamDataset::SetSpatialRef( poSRS ); + if( poSRS ) + m_oSRS = *poSRS; + else + m_oSRS.Clear(); if( m_poExternalDS ) - m_poExternalDS->SetProjection(pszProjection); + m_poExternalDS->SetSpatialRef(poSRS); InvalidateLabel(); return CE_None; } @@ -1948,19 +1945,25 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo ) if ((EQUAL( map_proj_name, "Equirectangular" )) || (EQUAL( map_proj_name, "SimpleCylindrical" )) ) { - oSRS.OGRSpatialReference::SetEquirectangular2 ( 0.0, center_lon, center_lat, 0, 0 ); + oSRS.SetEquirectangular2 ( 0.0, center_lon, center_lat, 0, 0 ); } else if (EQUAL( map_proj_name, "Orthographic" )) { - oSRS.OGRSpatialReference::SetOrthographic ( center_lat, center_lon, 0, 0 ); + oSRS.SetOrthographic ( center_lat, center_lon, 0, 0 ); } else if (EQUAL( map_proj_name, "Sinusoidal" )) { - oSRS.OGRSpatialReference::SetSinusoidal ( center_lon, 0, 0 ); + oSRS.SetSinusoidal ( center_lon, 0, 0 ); } else if (EQUAL( map_proj_name, "Mercator" )) { - oSRS.OGRSpatialReference::SetMercator ( center_lat, center_lon, scaleFactor, 0, 0 ); + oSRS.SetMercator ( center_lat, center_lon, scaleFactor, 0, 0 ); } else if (EQUAL( map_proj_name, "PolarStereographic" )) { - oSRS.OGRSpatialReference::SetPS ( center_lat, center_lon, scaleFactor, 0, 0 ); + oSRS.SetPS ( center_lat, center_lon, scaleFactor, 0, 0 ); } else if (EQUAL( map_proj_name, "TransverseMercator" )) { - oSRS.OGRSpatialReference::SetTM ( center_lat, center_lon, scaleFactor, 0, 0 ); + oSRS.SetTM ( center_lat, center_lon, scaleFactor, 0, 0 ); } else if (EQUAL( map_proj_name, "LambertConformal" )) { - oSRS.OGRSpatialReference::SetLCC ( first_std_parallel, second_std_parallel, center_lat, center_lon, 0, 0 ); + 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 { CPLDebug( "ISIS3", "Dataset projection %s is not supported. Continuing...", @@ -2015,7 +2018,8 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo ) else if ( (EQUAL( map_proj_name, "SimpleCylindrical" )) || (EQUAL( map_proj_name, "Orthographic" )) || (EQUAL( map_proj_name, "Stereographic" )) || - (EQUAL( map_proj_name, "Sinusoidal" )) ) { + (EQUAL( map_proj_name, "Sinusoidal" )) || + (EQUAL( map_proj_name, "PointPerspective" )) ) { // ISIS uses the spherical equation for these projections // so force a sphere. oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName, @@ -2053,10 +2057,8 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo ) } // translate back into a projection string. - char *pszResult = nullptr; - oSRS.exportToWkt( &pszResult ); - poDS->m_osProjection = pszResult; - CPLFree( pszResult ); + poDS->m_oSRS = oSRS; + poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); } /* END ISIS3 Label Read */ @@ -2480,10 +2482,8 @@ GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo ) if( oSRS2.importFromESRI( papszLines ) == OGRERR_NONE ) { poDS->m_aosAdditionalFiles.AddString( pszPrjFile ); - char *pszResult = nullptr; - oSRS2.exportToWkt( &pszResult ); - poDS->m_osProjection = pszResult; - CPLFree( pszResult ); + poDS->m_oSRS = oSRS2; + poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); } CSLDestroy( papszLines ); @@ -2670,7 +2670,7 @@ void ISIS3Dataset::BuildLabel() oPixels.Set( "Base", GetRasterBand(1)->GetOffset() ); oPixels.Set( "Multiplier", GetRasterBand(1)->GetScale() ); - OGRSpatialReference oSRS; + const OGRSpatialReference& oSRS = m_oSRS; if( !m_bUseSrcMapping ) { @@ -2688,11 +2688,10 @@ void ISIS3Dataset::BuildLabel() if( !m_osLongitudeDirection.empty() ) oMapping.Set( "LongitudeDirection", m_osLongitudeDirection ); } - else if( !m_bUseSrcMapping && !m_osProjection.empty() ) + else if( !m_bUseSrcMapping && !m_oSRS.IsEmpty() ) { oMapping.Add( "_type", "group" ); - oSRS.SetFromUserInput(m_osProjection); if( oSRS.IsProjected() || oSRS.IsGeographic() ) { const char* pszDatum = oSRS.GetAttrValue("DATUM"); @@ -2903,6 +2902,18 @@ void ISIS3Dataset::BuildLabel() oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0) ); } + else if( EQUAL(pszProjection, "Vertical Perspective") ) // PROJ 7 required + { + oMapping.Add( "ProjectionName", "PointPerspective" ); + oMapping.Add( "CenterLongitude", FixLong( + oSRS.GetNormProjParm("Longitude of topocentric origin", 0.0)) ); + oMapping.Add( "CenterLatitude", + oSRS.GetNormProjParm("Latitude of topocentric origin", 0.0) ); + // ISIS3 value is the distance from center of ellipsoid, in km + oMapping.Add( "Distance", + (oSRS.GetNormProjParm("Viewpoint height", 0.0) + oSRS.GetSemiMajor()) / 1000.0 ); + } + else { CPLError(CE_Warning, CPLE_NotSupported, @@ -2939,7 +2950,7 @@ void ISIS3Dataset::BuildLabel() oMapping.Add( "_type", "group" ); const double dfDegToMeter = oSRS.GetSemiMajor() * M_PI / 180.0; - if( !m_osProjection.empty() && oSRS.IsProjected() ) + if( !m_oSRS.IsEmpty() && oSRS.IsProjected() ) { const double dfLinearUnits = oSRS.GetLinearUnits(); // Maybe we should deal differently with non meter units ? @@ -2952,7 +2963,7 @@ void ISIS3Dataset::BuildLabel() oMapping.Add( "Scale/value", dfScale ); oMapping.Add( "Scale/unit", "pixels/degree" ); } - else if( !m_osProjection.empty() && oSRS.IsGeographic() ) + else if( !m_oSRS.IsEmpty() && oSRS.IsGeographic() ) { const double dfScale = 1.0 / m_adfGeoTransform[1]; const double dfRes = m_adfGeoTransform[1] * dfDegToMeter; @@ -4207,10 +4218,10 @@ GDALDataset* ISIS3Dataset::CreateCopy( const char *pszFilename, poDS->SetGeoTransform( adfGeoTransform ); } - if( poSrcDS->GetProjectionRef() != nullptr - && strlen(poSrcDS->GetProjectionRef()) > 0 ) + auto poSrcSRS = poSrcDS->GetSpatialRef(); + if( poSrcSRS ) { - poDS->SetProjection( poSrcDS->GetProjectionRef() ); + poDS->SetSpatialRef( poSrcSRS ); } for(int i=1;i<=nBands;i++)