diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e1993e125..ae9e9988bd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -56,6 +56,7 @@ release. override the timestamp style naming convention of the output cube with their own name; if not specified retains existing behavior [#5125](https://github.com/USGS-Astrogeology/ISIS3/issues/5162) - Added new parameters ONERROR, ERRORLOG, and ERRORLIST to mosrange to provide better control over error behavior and provide diagnostics when problems are encountered processing the input file list.[#3606](https://github.com/DOI-USGS/ISIS3/issues/3606) +- Adds PROJ into ISIS, and exposes the capability with a new class called IProj. [#5317](https://github.com/DOI-USGS/ISIS3/pull/5317) ### Deprecated diff --git a/isis/CMakeLists.txt b/isis/CMakeLists.txt index 0bb2e59308..806b98ff26 100644 --- a/isis/CMakeLists.txt +++ b/isis/CMakeLists.txt @@ -280,6 +280,7 @@ find_package(NN REQUIRED) find_package(OpenCV 3.1.0 REQUIRED) find_package(PCL REQUIRED) find_package(Protobuf 2.6.1 REQUIRED) +find_package(PROJ REQUIRED) find_package(Qwt 6 REQUIRED) find_package(SuperLU 4.3 REQUIRED) find_package(TIFF 4.0.0 REQUIRED) diff --git a/isis/src/base/apps/cam2map/cam2map.cpp b/isis/src/base/apps/cam2map/cam2map.cpp index 9f30e32896..217b601c29 100644 --- a/isis/src/base/apps/cam2map/cam2map.cpp +++ b/isis/src/base/apps/cam2map/cam2map.cpp @@ -30,7 +30,15 @@ namespace Isis { // Get the map projection file provided by the user Pvl userMap; - userMap.read(ui.GetFileName("MAP")); + if (ui.GetBoolean("USEPROJ")) { + PvlGroup mappingGroup("Mapping"); + mappingGroup.addKeyword(PvlKeyword("ProjectionName", "IProj")); + mappingGroup.addKeyword(PvlKeyword("ProjStr", ui.GetAsString("PROJString"))); + userMap.addGroup(mappingGroup); + } + else { + userMap.read(ui.GetFileName("MAP")); + } PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse); cam2map(&icube, userMap, userGrp, ui, log); diff --git a/isis/src/base/apps/cam2map/cam2map.xml b/isis/src/base/apps/cam2map/cam2map.xml index 567ab61422..3687b714c0 100644 --- a/isis/src/base/apps/cam2map/cam2map.xml +++ b/isis/src/base/apps/cam2map/cam2map.xml @@ -381,6 +381,28 @@ LONSEAM + + + boolean + false + + MAP + + + PROJString + + + + + + + string + A PROJ4 string that defines the projection + + This parameters is used instead of the map file to define the projection to + transform into. For help with PROJ strings visit https://proj.org/en/9.3/index.html + + diff --git a/isis/src/base/objs/IProj/IProj.cpp b/isis/src/base/objs/IProj/IProj.cpp new file mode 100644 index 0000000000..b6675f9d59 --- /dev/null +++ b/isis/src/base/objs/IProj/IProj.cpp @@ -0,0 +1,221 @@ +/** This is free and unencumbered software released into the public domain. +The authors of ISIS do not claim copyright on the contents of this file. +For more details about the LICENSE terms and the AUTHORS, you will +find files of those names at the top level of this repository. **/ + +/* SPDX-License-Identifier: CC0-1.0 */ +#include "IProj.h" + +#include + +#include "IException.h" +#include "TProjection.h" +#include "Pvl.h" +#include "PvlGroup.h" +#include "PvlKeyword.h" + +using namespace std; + +namespace Isis { + IProj::IProj(Pvl &label, bool allowDefaults) : + TProjection::TProjection(label) { + PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse); + if (!mapGroup.hasKeyword("ProjStr")) { + QString message = "No ProjStr keyword in mapping group, either add a ProjStr or select a different projection method"; + throw IException(IException::User, message, _FILEINFO_); + } + m_C = proj_context_create(); + + m_userOutputProjStr = new QString(mapGroup["ProjStr"]); + + /* Create a projection. */ + std::string projString = m_userOutputProjStr->toStdString(); + m_outputProj = proj_create(m_C, projString.c_str()); + + if (0 == m_outputProj) { + QString msg = "Unable to create projection from [" + QString(projString.c_str()) + "]"; + throw IException(IException::User, msg, _FILEINFO_); + } + + PJ *outputEllipsoid = proj_get_ellipsoid(m_C, m_outputProj); + + if (outputEllipsoid == 0) { + QString msg = "Unable to get ellipsoid from [" + *m_userOutputProjStr + "]. " + + "Please add a radii definition to your proj string."; + throw IException(IException::User, msg, _FILEINFO_); + } + + int res = proj_ellipsoid_get_parameters(m_C, outputEllipsoid, + &m_equatorialRadius, + &m_polarRadius, + nullptr, + nullptr); + proj_destroy(outputEllipsoid); + + if (res == 0) { + QString msg = "Unable to get ellipsoid information from [" + *m_userOutputProjStr + "]"; + throw IException(IException::User, msg, _FILEINFO_); + } + + m_llaProj = proj_crs_get_geodetic_crs(m_C, m_outputProj); + + if (0 == m_llaProj) { + QString msg = "Unable to create projection from []"; + throw IException(IException::Programmer, msg, _FILEINFO_); + } + + m_llaProj2outputProj = proj_create_crs_to_crs_from_pj(m_C, m_llaProj, m_outputProj, 0, 0); + if (0 == m_llaProj2outputProj) { + QString msg = "Unable to create transform from [" + QString(projString.c_str()) + "]"; + throw IException(IException::User, msg, _FILEINFO_); + } + } + + //! Destroys the IProj object + IProj::~IProj() { + delete m_userOutputProjStr; + proj_destroy(m_llaProj); + proj_destroy(m_outputProj); + proj_destroy(m_llaProj2outputProj); + proj_context_destroy(m_C); + } + + /** + * Returns the name of the map projection, "Proj" + * + * @return QString Name of projection, "Proj" + */ + QString IProj::Name() const { + return "Proj"; + } + + /** + * This function returns the keywords that this projection uses. + * + * We also include the generated PROJ string + * + * @return PvlGroup The keywords that this projection uses + */ + PvlGroup IProj::Mapping() { + PvlGroup mapping = TProjection::Mapping(); + + mapping.addKeyword(PvlKeyword("EquatorialRadius", toString(m_equatorialRadius, 15), "meters"), PvlContainer::InsertMode::Replace); + mapping.addKeyword(PvlKeyword("PolarRadius", toString(m_polarRadius, 15), "meters"), PvlContainer::InsertMode::Replace); + mapping.addKeyword(PvlKeyword("LatitudeType", "Planetographic"), PvlContainer::InsertMode::Replace); + mapping.addKeyword(PvlKeyword("LongitudeDomain", "180"), PvlContainer::InsertMode::Replace); + + mapping += PvlKeyword("ProjStr", *m_userOutputProjStr); + + return mapping; + } + + /** + * Returns the version of the map projection + * + * + * @return QString Version number + */ + QString IProj::Version() const { + return "1.0"; + } + + bool IProj::SetGround(const double lat, const double lon) { + m_longitude = lon; + m_latitude = lat; + + PJ_COORD c_in; + c_in.lpz.lam = m_longitude; + // Convert to ographic as proj defaults to operating in ographic latitudes + c_in.lpz.phi = ToPlanetographic(m_latitude); + + PJ_COORD c_out = proj_trans(m_llaProj2outputProj, PJ_FWD, c_in); + SetComputedXY(c_out.xy.x, c_out.xy.y); + m_good = true; + return m_good; + } + + bool IProj::SetCoordinate(const double x, const double y) { + SetXY(x, y); + + PJ_COORD c_in; + c_in.xy.x = x; + c_in.xy.y = y; + + PJ_COORD c_out = proj_trans(m_llaProj2outputProj, PJ_INV, c_in); + + m_longitude = c_out.lpz.lam; + // Convert back to ocentric as ISIS defaults to operating in ocentric latitudes + m_latitude = ToPlanetocentric(c_out.lpz.phi); + m_good = true; + return m_good; + } + + /** + * This method is used to determine the x/y range which completely covers the + * area of interest specified by the lat/lon range. The latitude/longitude + * range may be obtained from the labels. The purpose of this method is to + * return the x/y range so it can be used to compute how large a map may need + * to be. For example, how big a piece of paper is needed or how large of an + * image needs to be created. The method may fail as indicated by its return + * value. This currently mimics the sinusoidal projection's XYRange check + * and should be made more robust for use in proj. This will likely be a + * method that walks the boundary of the projection + * + * @param minX Minimum x projection coordinate which covers the latitude + * longitude range specified in the labels. + * + * @param maxX Maximum x projection coordinate which covers the latitude + * longitude range specified in the labels. + * + * @param minY Minimum y projection coordinate which covers the latitude + * longitude range specified in the labels. + * + * @param maxY Maximum y projection coordinate which covers the latitude + * longitude range specified in the labels. + * + * @return bool Indicates whether the method was successful. + */ + bool IProj::XYRange(double &minX, double &maxX, + double &minY, double &maxY) { + // Check the corners of the lat/lon range + XYRangeCheck(m_minimumLatitude, m_minimumLongitude); + XYRangeCheck(m_maximumLatitude, m_minimumLongitude); + XYRangeCheck(m_minimumLatitude, m_maximumLongitude); + XYRangeCheck(m_maximumLatitude, m_maximumLongitude); + + // If the latitude crosses the equator check there + if ((m_minimumLatitude < 0.0) && (m_maximumLatitude > 0.0)) { + XYRangeCheck(0.0, m_minimumLongitude); + XYRangeCheck(0.0, m_maximumLongitude); + } + + // Make sure everything is ordered + if (m_minimumX >= m_maximumX) return false; + if (m_minimumY >= m_maximumY) return false; + + // Return X/Y min/maxs + minX = m_minimumX; + maxX = m_maximumX; + minY = m_minimumY; + maxY = m_maximumY; + return true; + } + + /** + * This is the function that is called in order to instantiate a + * Sinusoidal object. + * + * @param lab Cube labels with appropriate Mapping information. + * + * @param allowDefaults Indicates whether CenterLongitude are allowed to + * be computed using the middle of the longitude + * range specified in the labels. + * + * @return @b Isis::Projection* Pointer to a Sinusoidal projection object. + */ + extern "C" Isis::TProjection *IProjPlugin(Isis::Pvl &lab, + bool allowDefaults) + { + return new Isis::IProj(lab, allowDefaults); + } +} \ No newline at end of file diff --git a/isis/src/base/objs/IProj/IProj.h b/isis/src/base/objs/IProj/IProj.h new file mode 100644 index 0000000000..0c46a00e70 --- /dev/null +++ b/isis/src/base/objs/IProj/IProj.h @@ -0,0 +1,60 @@ +#ifndef IProj_h +#define IProj_h +/** This is free and unencumbered software released into the public domain. +The authors of ISIS do not claim copyright on the contents of this file. +For more details about the LICENSE terms and the AUTHORS, you will +find files of those names at the top level of this repository. **/ + +/* SPDX-License-Identifier: CC0-1.0 */ +#include + +#include "TProjection.h" + +namespace Isis { + class Pvl; + class PvlGroup; + /** + * @brief Proj Map Projection + * + * This class provides methods for the forward and inverse projection for + * any ISIS map file through PROJ. We take a map file and convert it into + * a PROJ string, that string is then fed into the PROJ projection engine. + * + * + * Please see the TProjection/Projection class for a full accounting of all the methods + * available. + * + * @ingroup MapProjection + * + * @author 2023-10-17 Adam Paquette + * + * @internal + * @history 2023-10-17 Adam Paquette - Built Class + */ + + class IProj : public TProjection { + public: + IProj(Pvl &label, bool allowDefaults = false); + ~IProj(); + + QString Name() const; + QString Version() const; + + PvlGroup Mapping(); + + bool SetGround(const double lat, const double lon); + bool SetCoordinate(const double x, const double y); + + bool XYRange(double &minX, double &maxX, + double &minY, double &maxY); + + private: + QString *m_userOutputProjStr; + PJ_CONTEXT *m_C; + PJ *m_llaProj; + PJ *m_outputProj; + PJ *m_llaProj2outputProj; + }; +}; + +#endif diff --git a/isis/src/base/objs/IProj/Projection.plugin b/isis/src/base/objs/IProj/Projection.plugin new file mode 100644 index 0000000000..ac0d6ed06b --- /dev/null +++ b/isis/src/base/objs/IProj/Projection.plugin @@ -0,0 +1,5 @@ +Group = IProj + Library = IProj + Routine = IProjPlugin + Keyword = ProjStr +EndGroup diff --git a/isis/src/qisis/objs/AdvancedTrackTool/AdvancedTrackTool.cpp b/isis/src/qisis/objs/AdvancedTrackTool/AdvancedTrackTool.cpp index 3251048cd8..fa58f7ed90 100644 --- a/isis/src/qisis/objs/AdvancedTrackTool/AdvancedTrackTool.cpp +++ b/isis/src/qisis/objs/AdvancedTrackTool/AdvancedTrackTool.cpp @@ -581,7 +581,7 @@ namespace Isis { setText(QString::number(TProjection::To180Domain(lon), 'f', 15)); p_tableWin->table()->item(row, getIndex("360 Positive West Longitude"))-> setText(QString::number(wlon, 'f', 15)); - p_tableWin->table()->item(row, getIndex("180 Positive East Longitude"))-> + p_tableWin->table()->item(row, getIndex("180 Positive West Longitude"))-> setText(QString::number(TProjection::To180Domain(wlon), 'f', 15)); p_tableWin->table()->item(row, getIndex("Local Radius"))->setText(QString::number(radius, 'f', 15)); }