Skip to content

Commit

Permalink
ISIS3: add support for PointPerspective projection (refs OSGeo#1856)
Browse files Browse the repository at this point in the history
  • Loading branch information
@rouault
rouault committed Oct 2, 2019
1 parent 64d6ffe commit e6f4c66
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 43 deletions.
43 changes: 43 additions & 0 deletions autotest/gdrivers/data/isis3_pointperspective.cub
View file
Original file line number Diff line number Diff line change
@@ -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
29 changes: 29 additions & 0 deletions autotest/gdrivers/isis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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')
97 changes: 54 additions & 43 deletions gdal/frmts/pds/isis3dataset.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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...",
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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 );
Expand Down Expand Up @@ -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 )
{
Expand All @@ -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");
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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 ?
Expand All @@ -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;
Expand Down Expand Up @@ -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++)
Expand Down

0 comments on commit e6f4c66

Please sign in to comment.