Skip to content

Commit

Permalink
Merge pull request #2847 from rouault/spherical_ocentric
Browse files Browse the repository at this point in the history 
Add support for GeodeticCRS using a Spherical ocentric coordinate system
  • Loading branch information
@rouault
rouault authored Sep 14, 2021
2 parents 884cff4 + 078952e commit 92ca1a9
Show file tree
Hide file tree
Showing 17 changed files with 844 additions and 185 deletions.
4 changes: 4 additions & 0 deletions include/proj/coordinateoperation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1385,6 +1385,10 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
createGeographicGeocentric(const crs::CRSNNPtr &sourceCRS,
const crs::CRSNNPtr &targetCRS);

PROJ_INTERNAL static ConversionNNPtr
createGeographicGeocentricLatitude(const crs::CRSNNPtr &sourceCRS,
const crs::CRSNNPtr &targetCRS);

//! @endcond

protected:
Expand Down
5 changes: 5 additions & 0 deletions include/proj/coordinatesystem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,11 @@ class PROJ_GCC_DLL SphericalCS final : public CoordinateSystem {
const CoordinateSystemAxisNNPtr &axis2,
const CoordinateSystemAxisNNPtr &axis3);

PROJ_DLL static SphericalCSNNPtr
create(const util::PropertyMap &properties,
const CoordinateSystemAxisNNPtr &axis1,
const CoordinateSystemAxisNNPtr &axis2);

protected:
PROJ_INTERNAL explicit SphericalCS(
const std::vector<CoordinateSystemAxisNNPtr> &axisIn);
Expand Down
13 changes: 8 additions & 5 deletions include/proj/crs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,6 +270,8 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS,

PROJ_DLL bool isGeocentric() PROJ_PURE_DECL;

PROJ_DLL bool isSphericalPlanetocentric() PROJ_PURE_DECL;

PROJ_DLL static GeodeticCRSNNPtr
create(const util::PropertyMap &properties,
const datum::GeodeticReferenceFrameNNPtr &datum,
Expand Down Expand Up @@ -308,6 +310,9 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS,
PROJ_INTERNAL void addGeocentricUnitConversionIntoPROJString(
io::PROJStringFormatter *formatter) const;

PROJ_INTERNAL void
addAngularUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter) const;

PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter)
const override; // throw(io::FormattingException)

Expand Down Expand Up @@ -400,12 +405,10 @@ class PROJ_GCC_DLL GeographicCRS : public GeodeticCRS {

PROJ_PRIVATE :
//! @cond Doxygen_Suppress
PROJ_INTERNAL void
addAngularUnitConvertAndAxisSwap(
io::PROJStringFormatter *formatter) const;

PROJ_INTERNAL void _exportToPROJString(io::PROJStringFormatter *formatter)
const override; // throw(FormattingException)
PROJ_INTERNAL void
_exportToPROJString(io::PROJStringFormatter *formatter)
const override; // throw(FormattingException)

PROJ_INTERNAL void _exportToJSON(io::JSONFormatter *formatter)
const override; // throw(FormattingException)
Expand Down
2 changes: 2 additions & 0 deletions scripts/reference_exported_symbols.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@ osgeo::proj::crs::GeodeticCRS::ellipsoid() const
osgeo::proj::crs::GeodeticCRS::~GeodeticCRS()
osgeo::proj::crs::GeodeticCRS::identify(std::shared_ptr<osgeo::proj::io::AuthorityFactory> const&) const
osgeo::proj::crs::GeodeticCRS::isGeocentric() const
osgeo::proj::crs::GeodeticCRS::isSphericalPlanetocentric() const
osgeo::proj::crs::GeodeticCRS::primeMeridian() const
osgeo::proj::crs::GeodeticCRS::velocityModel() const
osgeo::proj::crs::GeographicCRS::coordinateSystem() const
Expand Down Expand Up @@ -225,6 +226,7 @@ osgeo::proj::cs::OrdinalCS::create(osgeo::proj::util::PropertyMap const&, std::v
osgeo::proj::cs::OrdinalCS::~OrdinalCS()
osgeo::proj::cs::ParametricCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
osgeo::proj::cs::ParametricCS::~ParametricCS()
osgeo::proj::cs::SphericalCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
osgeo::proj::cs::SphericalCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
osgeo::proj::cs::SphericalCS::~SphericalCS()
osgeo::proj::cs::TemporalCountCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn<std::shared_ptr<osgeo::proj::cs::CoordinateSystemAxis> > const&)
Expand Down
33 changes: 18 additions & 15 deletions src/apps/cs2cs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,10 +59,10 @@

static PJ *transformation = nullptr;

static bool srcIsGeog = false;
static bool srcIsLongLat = false;
static double srcToRadians = 0.0;

static bool destIsGeog = false;
static bool destIsLongLat = false;
static double destToRadians = 0.0;
static bool destIsLatLong = false;

Expand Down Expand Up @@ -158,7 +158,7 @@ static void process(FILE *fid)

if (data.u != HUGE_VAL) {

if (srcIsGeog) {
if (srcIsLongLat) {
/* dmstor gives values to radians. Convert now to the SRS unit
*/
data.u /= srcToRadians;
Expand All @@ -179,7 +179,7 @@ static void process(FILE *fid)
if (data.u == HUGE_VAL) /* error output */
fputs(oterr, stdout);

else if (destIsGeog && !oform) { /*ascii DMS output */
else if (destIsLongLat && !oform) { /*ascii DMS output */

// rtodms() expect radians: convert from the output SRS unit
data.u *= destToRadians;
Expand All @@ -206,7 +206,7 @@ static void process(FILE *fid)
}

} else { /* x-y or decimal degree ascii output */
if (destIsGeog) {
if (destIsLongLat) {
data.v *= destToRadians * RAD_TO_DEG;
data.u *= destToRadians * RAD_TO_DEG;
}
Expand Down Expand Up @@ -239,15 +239,15 @@ static void process(FILE *fid)
/************************************************************************/

static PJ *instantiate_crs(const std::string &definition,
bool &isGeog, double &toRadians,
bool &isLongLatCS, double &toRadians,
bool &isLatFirst) {
PJ *crs = proj_create(nullptr,
pj_add_type_crs_if_needed(definition).c_str());
if (!crs) {
return nullptr;
}

isGeog = false;
isLongLatCS = false;
toRadians = 0.0;
isLatFirst = false;

Expand All @@ -259,11 +259,11 @@ static PJ *instantiate_crs(const std::string &definition,
type = proj_get_type(crs);
}
if (type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
type == PJ_TYPE_GEOGRAPHIC_3D_CRS) {
type == PJ_TYPE_GEOGRAPHIC_3D_CRS ||
type == PJ_TYPE_GEODETIC_CRS) {
auto cs = proj_crs_get_coordinate_system(nullptr, crs);
assert(cs);

isGeog = true;
const char *axisName = "";
proj_cs_get_axis_info(nullptr, cs, 0,
&axisName, // name,
Expand All @@ -277,6 +277,9 @@ static PJ *instantiate_crs(const std::string &definition,
isLatFirst =
NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
std::string::npos;
isLongLatCS = isLatFirst ||
NS_PROJ::internal::ci_find(std::string(axisName), "longitude") !=
std::string::npos;

proj_destroy(cs);
}
Expand Down Expand Up @@ -736,7 +739,7 @@ int main(int argc, char **argv) {
PJ *src = nullptr;
if (!fromStr.empty()) {
bool ignored;
src = instantiate_crs(fromStr, srcIsGeog,
src = instantiate_crs(fromStr, srcIsLongLat,
srcToRadians, ignored);
if (!src) {
emess(3, "cannot instantiate source coordinate system");
Expand All @@ -745,7 +748,7 @@ int main(int argc, char **argv) {

PJ *dst = nullptr;
if (!toStr.empty()) {
dst = instantiate_crs(toStr, destIsGeog,
dst = instantiate_crs(toStr, destIsLongLat,
destToRadians, destIsLatLong);
if (!dst) {
emess(3, "cannot instantiate target coordinate system");
Expand All @@ -760,7 +763,7 @@ int main(int argc, char **argv) {
emess(3,
"missing target CRS and source CRS is not a projected CRS");
}
destIsGeog = true;
destIsLongLat = true;
} else if (fromStr.empty()) {
assert(dst);
bool ignored;
Expand All @@ -770,7 +773,7 @@ int main(int argc, char **argv) {
emess(3,
"missing source CRS and target CRS is not a projected CRS");
}
srcIsGeog = true;
srcIsLongLat = true;
}
proj_destroy(src);
proj_destroy(dst);
Expand Down Expand Up @@ -835,13 +838,13 @@ int main(int argc, char **argv) {
}

/* set input formatting control */
if (!srcIsGeog)
if (!srcIsLongLat)
informat = strtod;
else {
informat = dmstor;
}

if (!destIsGeog && !oform)
if (!destIsLongLat && !oform)
oform = "%.2f";

/* process input file list */
Expand Down
21 changes: 21 additions & 0 deletions src/iso19111/coordinatesystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -667,6 +667,27 @@ SphericalCSNNPtr SphericalCS::create(const util::PropertyMap &properties,

// ---------------------------------------------------------------------------

/** \brief Instantiate a SphericalCS with 2 axis.
*
* This is an extension to ISO19111 to support (planet)-ocentric CS with
* geocentric latitude.
*
* @param properties See \ref general_properties.
* @param axis1 The first axis.
* @param axis2 The second axis.
* @return a new SphericalCS.
*/
SphericalCSNNPtr SphericalCS::create(const util::PropertyMap &properties,
const CoordinateSystemAxisNNPtr &axis1,
const CoordinateSystemAxisNNPtr &axis2) {
std::vector<CoordinateSystemAxisNNPtr> axis{axis1, axis2};
auto cs(SphericalCS::nn_make_shared<SphericalCS>(axis));
cs->setProperties(properties);
return cs;
}

// ---------------------------------------------------------------------------

//! @cond Doxygen_Suppress
EllipsoidalCS::~EllipsoidalCS() = default;
//! @endcond
Expand Down
Loading

0 comments on commit 92ca1a9

Please sign in to comment.