From 9150723471a128bdd24341deb0df82115c68ba5f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 2 Oct 2019 22:57:09 +0200 Subject: [PATCH 1/4] Add support for Vertical Perspective projection (refs #1856) --- autotest/osr/osr_basic.py | 11 +++ gdal/ogr/ogr_spatialref.h | 9 ++ gdal/ogr/ogr_srs_api.h | 9 ++ gdal/ogr/ogrspatialreference.cpp | 85 ++++++++++++++++- gdal/swig/include/osr.i | 11 +++ gdal/swig/python/extensions/osr_wrap.cpp | 114 +++++++++++++++++++++++ gdal/swig/python/osgeo/osr.py | 5 + 7 files changed, 239 insertions(+), 5 deletions(-) diff --git a/autotest/osr/osr_basic.py b/autotest/osr/osr_basic.py index 4a953dfdc8aa..9f7476cb5867 100755 --- a/autotest/osr/osr_basic.py +++ b/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 + diff --git a/gdal/ogr/ogr_spatialref.h b/gdal/ogr/ogr_spatialref.h index 73c397c6c2f0..a79b208ce129 100644 --- a/gdal/ogr/ogr_spatialref.h +++ b/gdal/ogr/ogr_spatialref.h @@ -638,6 +638,15 @@ class CPL_DLL OGRSpatialReference /** Spherical, Cross-track, Height */ OGRErr SetSCH( double dfPegLat, double dfPegLong, double dfPegHeading, double dfPegHgt); + + /** Vertical Perspective / Near-sided Perspective */ + OGRErr SetVerticalPerspective( double dfTopoOriginLat, + double dfTopoOriginLon, + double dfTopoOriginHeight, + double dfViewPointHeight, + double dfFalseEasting, + double dfFalseNorthing); + /** State Plane */ OGRErr SetStatePlane( int nZone, int bNAD83 = TRUE, const char *pszOverrideUnitName = nullptr, diff --git a/gdal/ogr/ogr_srs_api.h b/gdal/ogr/ogr_srs_api.h index 4c0d1413b65a..65a940c6f8bd 100644 --- a/gdal/ogr/ogr_srs_api.h +++ b/gdal/ogr/ogr_srs_api.h @@ -945,6 +945,15 @@ OGRErr CPL_DLL OSRSetSCH( OGRSpatialReferenceH hSRS, double dfPegLat, double dfPegLong, double dfPegHeading, double dfPegHgt); +/** Vertical Perspective / Near-sided Perspective */ +OGRErr CPL_DLL OSRSetVerticalPerspective( OGRSpatialReferenceH hSRS, + double dfTopoOriginLat, + double dfTopoOriginLon, + double dfTopoOriginHeight, + double dfViewPointHeight, + double dfFalseEasting, + double dfFalseNorthing); + double CPL_DLL OSRCalcInvFlattening( double dfSemiMajor, double dfSemiMinor ); double CPL_DLL OSRCalcSemiMinorFromInvFlattening( double dfSemiMajor, double dfInvFlattening ); diff --git a/gdal/ogr/ogrspatialreference.cpp b/gdal/ogr/ogrspatialreference.cpp index 21b68eac95ad..9196d01f7b10 100644 --- a/gdal/ogr/ogrspatialreference.cpp +++ b/gdal/ogr/ogrspatialreference.cpp @@ -92,6 +92,7 @@ struct OGRSpatialReference::Private CPLString m_osAreaName{}; bool m_bNodesChanged = false; + bool m_bNodesWKT2 = false; OGR_SRSNode *m_poRoot = nullptr; double dfFromGreenwich = 0.0; @@ -274,9 +275,22 @@ void OGRSpatialReference::Private::refreshRootFromProjObj() } aosOptions.SetNameValue("STRICT", "NO"); - const char* pszWKT = proj_as_wkt(getPROJContext(), - m_pj_crs, m_bMorphToESRI ? PJ_WKT1_ESRI : PJ_WKT1_GDAL, - aosOptions.List()); + const char* pszWKT; + { + CPLErrorStateBackuper oErrorStateBackuper; + CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler); + pszWKT = proj_as_wkt(getPROJContext(), + m_pj_crs, m_bMorphToESRI ? PJ_WKT1_ESRI : PJ_WKT1_GDAL, + aosOptions.List()); + m_bNodesWKT2 = false; + } + if( !m_bMorphToESRI && pszWKT == nullptr ) + { + pszWKT = proj_as_wkt(getPROJContext(), m_pj_crs, PJ_WKT2_2018, + aosOptions.List()); + m_bNodesWKT2 = true; + } + CPLPopErrorHandler(); if( pszWKT ) { auto root = new OGR_SRSNode(); @@ -1077,7 +1091,13 @@ const char *OGRSpatialReference::GetAttrValue( const char * pszNodeName, { const OGR_SRSNode *poNode = GetAttrNode( pszNodeName ); if( poNode == nullptr ) + { + if( d->m_bNodesWKT2 && EQUAL(pszNodeName, "PROJECTION") ) + { + return GetAttrValue("METHOD", iAttr); + } return nullptr; + } if( iAttr < 0 || iAttr >= poNode->GetChildCount() ) return nullptr; @@ -4864,7 +4884,7 @@ int OGRSpatialReference::FindProjParm( const char *pszParameter, const OGR_SRSNode *poParameter = poPROJCS->GetChild(iChild); if( EQUAL(poParameter->GetValue(), "PARAMETER") - && poParameter->GetChildCount() == 2 + && poParameter->GetChildCount() >= 2 && EQUAL(poPROJCS->GetChild(iChild)->GetChild(0)->GetValue(), pszParameter) ) { @@ -4924,7 +4944,7 @@ double OGRSpatialReference::GetProjParm( const char * pszName, /* -------------------------------------------------------------------- */ /* Find the desired parameter. */ /* -------------------------------------------------------------------- */ - const OGR_SRSNode *poPROJCS = GetAttrNode( "PROJCS" ); + const OGR_SRSNode *poPROJCS = GetAttrNode( d->m_bNodesWKT2 ? "CONVERSION" : "PROJCS" ); if( poPROJCS == nullptr ) { if( pnErr != nullptr ) @@ -7260,6 +7280,61 @@ OGRErr OSRSetSCH( OGRSpatialReferenceH hSRS, dfPegLat, dfPegLong, dfPegHeading, dfPegHgt ); } + +/************************************************************************/ +/* SetVerticalPerspective() */ +/************************************************************************/ + +OGRErr OGRSpatialReference::SetVerticalPerspective( double dfTopoOriginLat, + double dfTopoOriginLon, + double dfTopoOriginHeight, + double dfViewPointHeight, + double dfFalseEasting, + double dfFalseNorthing ) +{ +#if PROJ_VERSION_MAJOR >= 7 + return d->replaceConversionAndUnref( + proj_create_conversion_vertical_perspective( + d->getPROJContext(), + dfTopoOriginLat, dfTopoOriginLon, + dfTopoOriginHeight, dfViewPointHeight, + dfFalseEasting, dfFalseNorthing, + nullptr, 0, nullptr, 0)); +#else + CPL_IGNORE_RET_VAL(dfTopoOriginHeight); // ignored by PROJ + + OGRSpatialReference oSRS; + CPLString oProj4String; + oProj4String.Printf( + "+proj=nsper +lat_0=%.18g +lon_0=%.18g +h=%.18g +x_0=%.18g +y_0=%.18g", + dfTopoOriginLat, dfTopoOriginLon, dfViewPointHeight, + dfFalseEasting, dfFalseNorthing); + oSRS.SetFromUserInput(oProj4String); + return d->replaceConversionAndUnref( + proj_crs_get_coordoperation(d->getPROJContext(), oSRS.d->m_pj_crs)); +#endif +} + +/************************************************************************/ +/* OSRSetVerticalPerspective() */ +/************************************************************************/ + +OGRErr OSRSetVerticalPerspective( OGRSpatialReferenceH hSRS, + double dfTopoOriginLat, + double dfTopoOriginLon, + double dfTopoOriginHeight, + double dfViewPointHeight, + double dfFalseEasting, + double dfFalseNorthing ) + +{ + VALIDATE_POINTER1( hSRS, "OSRSetVerticalPerspective", OGRERR_FAILURE ); + + return ToPointer(hSRS)->SetVerticalPerspective( + dfTopoOriginLat, dfTopoOriginLon, dfTopoOriginHeight, + dfViewPointHeight, dfFalseEasting, dfFalseNorthing ); +} + /************************************************************************/ /* SetAuthority() */ /************************************************************************/ diff --git a/gdal/swig/include/osr.i b/gdal/swig/include/osr.i index 4f4300f656b2..08b90ee386c8 100644 --- a/gdal/swig/include/osr.i +++ b/gdal/swig/include/osr.i @@ -885,6 +885,17 @@ public: return OSRSetVDG( self, clong, fe, fn ); } +%feature( "kwargs" ) SetVerticalPerspective; + OGRErr SetVerticalPerspective( double topoOriginLat, + double topoOriginLon, + double topoOriginHeight, + double viewPointHeight, + double fe, double fn ) + { + return OSRSetVerticalPerspective( self, + topoOriginLat, topoOriginLon, topoOriginHeight, viewPointHeight, fe, fn ); + } + OGRErr SetWellKnownGeogCS( const char *name ) { return OSRSetWellKnownGeogCS( self, name ); } diff --git a/gdal/swig/python/extensions/osr_wrap.cpp b/gdal/swig/python/extensions/osr_wrap.cpp index 5baf5e65cf57..7439b55b68bc 100644 --- a/gdal/swig/python/extensions/osr_wrap.cpp +++ b/gdal/swig/python/extensions/osr_wrap.cpp @@ -4228,6 +4228,10 @@ SWIGINTERN OGRErr OSRSpatialReferenceShadow_SetTMSO(OSRSpatialReferenceShadow *s SWIGINTERN OGRErr OSRSpatialReferenceShadow_SetVDG(OSRSpatialReferenceShadow *self,double clong,double fe,double fn){ return OSRSetVDG( self, clong, fe, fn ); } +SWIGINTERN OGRErr OSRSpatialReferenceShadow_SetVerticalPerspective(OSRSpatialReferenceShadow *self,double topoOriginLat,double topoOriginLon,double topoOriginHeight,double viewPointHeight,double fe,double fn){ + return OSRSetVerticalPerspective( self, + topoOriginLat, topoOriginLon, topoOriginHeight, viewPointHeight, fe, fn ); + } SWIGINTERN OGRErr OSRSpatialReferenceShadow_SetWellKnownGeogCS(OSRSpatialReferenceShadow *self,char const *name){ return OSRSetWellKnownGeogCS( self, name ); } @@ -11920,6 +11924,115 @@ SWIGINTERN PyObject *_wrap_SpatialReference_SetVDG(PyObject *SWIGUNUSEDPARM(self } +SWIGINTERN PyObject *_wrap_SpatialReference_SetVerticalPerspective(PyObject *SWIGUNUSEDPARM(self), PyObject *args, PyObject *kwargs) { + PyObject *resultobj = 0; int bLocalUseExceptionsCode = bUseExceptions; + OSRSpatialReferenceShadow *arg1 = (OSRSpatialReferenceShadow *) 0 ; + double arg2 ; + double arg3 ; + double arg4 ; + double arg5 ; + double arg6 ; + double arg7 ; + void *argp1 = 0 ; + int res1 = 0 ; + double val2 ; + int ecode2 = 0 ; + double val3 ; + int ecode3 = 0 ; + double val4 ; + int ecode4 = 0 ; + double val5 ; + int ecode5 = 0 ; + double val6 ; + int ecode6 = 0 ; + double val7 ; + int ecode7 = 0 ; + PyObject * obj0 = 0 ; + PyObject * obj1 = 0 ; + PyObject * obj2 = 0 ; + PyObject * obj3 = 0 ; + PyObject * obj4 = 0 ; + PyObject * obj5 = 0 ; + PyObject * obj6 = 0 ; + char * kwnames[] = { + (char *) "self",(char *) "topoOriginLat",(char *) "topoOriginLon",(char *) "topoOriginHeight",(char *) "viewPointHeight",(char *) "fe",(char *) "fn", NULL + }; + OGRErr result; + + if (!PyArg_ParseTupleAndKeywords(args,kwargs,(char *)"OOOOOOO:SpatialReference_SetVerticalPerspective",kwnames,&obj0,&obj1,&obj2,&obj3,&obj4,&obj5,&obj6)) SWIG_fail; + res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_OSRSpatialReferenceShadow, 0 | 0 ); + if (!SWIG_IsOK(res1)) { + SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "1"" of type '" "OSRSpatialReferenceShadow *""'"); + } + arg1 = reinterpret_cast< OSRSpatialReferenceShadow * >(argp1); + ecode2 = SWIG_AsVal_double(obj1, &val2); + if (!SWIG_IsOK(ecode2)) { + SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "2"" of type '" "double""'"); + } + arg2 = static_cast< double >(val2); + ecode3 = SWIG_AsVal_double(obj2, &val3); + if (!SWIG_IsOK(ecode3)) { + SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "3"" of type '" "double""'"); + } + arg3 = static_cast< double >(val3); + ecode4 = SWIG_AsVal_double(obj3, &val4); + if (!SWIG_IsOK(ecode4)) { + SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "4"" of type '" "double""'"); + } + arg4 = static_cast< double >(val4); + ecode5 = SWIG_AsVal_double(obj4, &val5); + if (!SWIG_IsOK(ecode5)) { + SWIG_exception_fail(SWIG_ArgError(ecode5), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "5"" of type '" "double""'"); + } + arg5 = static_cast< double >(val5); + ecode6 = SWIG_AsVal_double(obj5, &val6); + if (!SWIG_IsOK(ecode6)) { + SWIG_exception_fail(SWIG_ArgError(ecode6), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "6"" of type '" "double""'"); + } + arg6 = static_cast< double >(val6); + ecode7 = SWIG_AsVal_double(obj6, &val7); + if (!SWIG_IsOK(ecode7)) { + SWIG_exception_fail(SWIG_ArgError(ecode7), "in method '" "SpatialReference_SetVerticalPerspective" "', argument " "7"" of type '" "double""'"); + } + arg7 = static_cast< double >(val7); + { + if ( bUseExceptions ) { + ClearErrorState(); + } + result = (OGRErr)OSRSpatialReferenceShadow_SetVerticalPerspective(arg1,arg2,arg3,arg4,arg5,arg6,arg7); +#ifndef SED_HACKS + if ( bUseExceptions ) { + CPLErr eclass = CPLGetLastErrorType(); + if ( eclass == CE_Failure || eclass == CE_Fatal ) { + SWIG_exception( SWIG_RuntimeError, CPLGetLastErrorMsg() ); + } + } +#endif + } + { + /* %typemap(out) OGRErr */ + if ( result != 0 && bUseExceptions) { + const char* pszMessage = CPLGetLastErrorMsg(); + if( pszMessage[0] != '\0' ) + PyErr_SetString( PyExc_RuntimeError, pszMessage ); + else + PyErr_SetString( PyExc_RuntimeError, OGRErrMessages(result) ); + SWIG_fail; + } + } + { + /* %typemap(ret) OGRErr */ + if ( ReturnSame(resultobj == Py_None || resultobj == 0) ) { + resultobj = PyInt_FromLong( result ); + } + } + if ( ReturnSame(bLocalUseExceptionsCode) ) { CPLErr eclass = CPLGetLastErrorType(); if ( eclass == CE_Failure || eclass == CE_Fatal ) { Py_XDECREF(resultobj); SWIG_Error( SWIG_RuntimeError, CPLGetLastErrorMsg() ); return NULL; } } + return resultobj; +fail: + return NULL; +} + + SWIGINTERN PyObject *_wrap_SpatialReference_SetWellKnownGeogCS(PyObject *SWIGUNUSEDPARM(self), PyObject *args) { PyObject *resultobj = 0; int bLocalUseExceptionsCode = bUseExceptions; OSRSpatialReferenceShadow *arg1 = (OSRSpatialReferenceShadow *) 0 ; @@ -17321,6 +17434,7 @@ static PyMethodDef SwigMethods[] = { { (char *)"SpatialReference_SetTMG", (PyCFunction) _wrap_SpatialReference_SetTMG, METH_VARARGS | METH_KEYWORDS, (char *)"SpatialReference_SetTMG(SpatialReference self, double clat, double clong, double fe, double fn) -> OGRErr"}, { (char *)"SpatialReference_SetTMSO", (PyCFunction) _wrap_SpatialReference_SetTMSO, METH_VARARGS | METH_KEYWORDS, (char *)"SpatialReference_SetTMSO(SpatialReference self, double clat, double clong, double scale, double fe, double fn) -> OGRErr"}, { (char *)"SpatialReference_SetVDG", (PyCFunction) _wrap_SpatialReference_SetVDG, METH_VARARGS | METH_KEYWORDS, (char *)"SpatialReference_SetVDG(SpatialReference self, double clong, double fe, double fn) -> OGRErr"}, + { (char *)"SpatialReference_SetVerticalPerspective", (PyCFunction) _wrap_SpatialReference_SetVerticalPerspective, METH_VARARGS | METH_KEYWORDS, (char *)"SpatialReference_SetVerticalPerspective(SpatialReference self, double topoOriginLat, double topoOriginLon, double topoOriginHeight, double viewPointHeight, double fe, double fn) -> OGRErr"}, { (char *)"SpatialReference_SetWellKnownGeogCS", _wrap_SpatialReference_SetWellKnownGeogCS, METH_VARARGS, (char *)"SpatialReference_SetWellKnownGeogCS(SpatialReference self, char const * name) -> OGRErr"}, { (char *)"SpatialReference_SetFromUserInput", _wrap_SpatialReference_SetFromUserInput, METH_VARARGS, (char *)"SpatialReference_SetFromUserInput(SpatialReference self, char const * name) -> OGRErr"}, { (char *)"SpatialReference_CopyGeogCSFrom", _wrap_SpatialReference_CopyGeogCSFrom, METH_VARARGS, (char *)"SpatialReference_CopyGeogCSFrom(SpatialReference self, SpatialReference rhs) -> OGRErr"}, diff --git a/gdal/swig/python/osgeo/osr.py b/gdal/swig/python/osgeo/osr.py index 1ecd15ba88cc..3e2d4aaffb9a 100644 --- a/gdal/swig/python/osgeo/osr.py +++ b/gdal/swig/python/osgeo/osr.py @@ -811,6 +811,11 @@ def SetVDG(self, *args, **kwargs): return _osr.SpatialReference_SetVDG(self, *args, **kwargs) + def SetVerticalPerspective(self, *args, **kwargs): + """SetVerticalPerspective(SpatialReference self, double topoOriginLat, double topoOriginLon, double topoOriginHeight, double viewPointHeight, double fe, double fn) -> OGRErr""" + return _osr.SpatialReference_SetVerticalPerspective(self, *args, **kwargs) + + def SetWellKnownGeogCS(self, *args): """SetWellKnownGeogCS(SpatialReference self, char const * name) -> OGRErr""" return _osr.SpatialReference_SetWellKnownGeogCS(self, *args) From 022246895ff4ab8555a6c5c8860a02f890cab302 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 2 Oct 2019 22:58:24 +0200 Subject: [PATCH 2/4] GTiff and PAM: allow serializing WKT2 for SRS using non-WKT1 compatible projections such as Vertical Perspective (refs #1856) --- autotest/gcore/tiff_srs.py | 28 ++++++++++++++++++++++++++++ gdal/frmts/gtiff/geotiff.cpp | 35 +++++++++++++++++++++++++++++------ gdal/gcore/gdalpamdataset.cpp | 12 +++++++++++- 3 files changed, 68 insertions(+), 7 deletions(-) diff --git a/autotest/gcore/tiff_srs.py b/autotest/gcore/tiff_srs.py index 87895b90ce04..237ce0d0ae2a 100755 --- a/autotest/gcore/tiff_srs.py +++ b/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') diff --git a/gdal/frmts/gtiff/geotiff.cpp b/gdal/frmts/gtiff/geotiff.cpp index 2c5482764d01..c47058f173c7 100644 --- a/gdal/frmts/gtiff/geotiff.cpp +++ b/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,6 +17570,7 @@ 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 ); @@ -17567,10 +17578,22 @@ GTiffDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 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()); diff --git a/gdal/gcore/gdalpamdataset.cpp b/gdal/gcore/gdalpamdataset.cpp index 39bc4912df84..c2b78c06fb30 100644 --- a/gdal/gcore/gdalpamdataset.cpp +++ b/gdal/gcore/gdalpamdataset.cpp @@ -195,7 +195,17 @@ CPLXMLNode *GDALPamDataset::SerializeToXML( const char *pszUnused ) if( psPam->poSRS && !psPam->poSRS->IsEmpty() ) { char* pszWKT = nullptr; - psPam->poSRS->exportToWkt(&pszWKT); + { + CPLErrorStateBackuper oErrorStateBackuper; + CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler); + if( psPam->poSRS->exportToWkt(&pszWKT) != OGRERR_NONE ) + { + CPLFree(pszWKT); + pszWKT = nullptr; + const char* const apszOptions[] = { "FORMAT=WKT2", nullptr }; + psPam->poSRS->exportToWkt(&pszWKT, apszOptions); + } + } CPLXMLNode* psSRSNode = CPLCreateXMLElementAndValue( psDSTree, "SRS", pszWKT ); CPLFree(pszWKT); const auto& mapping = psPam->poSRS->GetDataAxisToSRSAxisMapping(); From 64d6ffeadd7b3ebda4aa1265ae821fa6ac290d7b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 2 Oct 2019 22:58:58 +0200 Subject: [PATCH 3/4] gdal_translate / gdalwarp / ogrct: allow dealing with non-WKT1 representable SRS (refs #1856) --- gdal/apps/gdal_translate_lib.cpp | 12 ++++- gdal/apps/gdalwarp_lib.cpp | 76 ++++++++++++++++++++------------ gdal/ogr/ogrct.cpp | 16 +++++-- 3 files changed, 72 insertions(+), 32 deletions(-) diff --git a/gdal/apps/gdal_translate_lib.cpp b/gdal/apps/gdal_translate_lib.cpp index dddd56f622c1..86f73839eeec 100644 --- a/gdal/apps/gdal_translate_lib.cpp +++ b/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 ); diff --git a/gdal/apps/gdalwarp_lib.cpp b/gdal/apps/gdalwarp_lib.cpp index 604754586b55..ed1132507089 100644 --- a/gdal/apps/gdalwarp_lib.cpp +++ b/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; diff --git a/gdal/ogr/ogrct.cpp b/gdal/ogr/ogrct.cpp index 7043ed86f6ca..28a4405ab18b 100644 --- a/gdal/ogr/ogrct.cpp +++ b/gdal/ogr/ogrct.cpp @@ -686,8 +686,13 @@ int OGRProjCT::Initialize( const OGRSpatialReference * poSourceIn, CPLDebug( "OGRCT", "Wrap at %g.", dfSourceWrapLong ); } - const char *pszCENTER_LONG = - poSRSSource ? poSRSSource->GetExtension( "GEOGCS", "CENTER_LONG" ) : nullptr; + const char *pszCENTER_LONG; + { + CPLErrorStateBackuper oErrorStateBackuper; + CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler); + pszCENTER_LONG = + poSRSSource ? poSRSSource->GetExtension( "GEOGCS", "CENTER_LONG" ) : nullptr; + } if( pszCENTER_LONG != nullptr ) { dfSourceWrapLong = CPLAtof(pszCENTER_LONG); @@ -701,7 +706,12 @@ int OGRProjCT::Initialize( const OGRSpatialReference * poSourceIn, CPLDebug( "OGRCT", "Wrap source at %g.", dfSourceWrapLong ); } - pszCENTER_LONG = poSRSTarget ? poSRSTarget->GetExtension( "GEOGCS", "CENTER_LONG" ) : nullptr; + { + CPLErrorStateBackuper oErrorStateBackuper; + CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler); + pszCENTER_LONG = poSRSTarget ? + poSRSTarget->GetExtension( "GEOGCS", "CENTER_LONG" ) : nullptr; + } if( pszCENTER_LONG != nullptr ) { dfTargetWrapLong = CPLAtof(pszCENTER_LONG); From e6f4c6697477ce23c3debabe21fae07d5c74c84c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 2 Oct 2019 22:59:28 +0200 Subject: [PATCH 4/4] 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++)