Skip to content

GTiff SRS Propagation #5736

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Kelvinrr merged 5 commits into from
Mar 7, 2025
Merged

GTiff SRS Propagation #5736

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
4a344dc
Set the SpatialReference for GTiffs that use the proj projection plugin
acpaquette Feb 27, 2025
ac35682
Added changelog entry
acpaquette Feb 27, 2025
a1cb7a9
Added test for opening a non-gdal created ISIS cube and writing a map…
acpaquette Mar 3, 2025
8da5ed9
Fixed opening/reopening with blobs
acpaquette Mar 4, 2025
29bc4a8
Fixed GTiff test name
acpaquette Mar 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ release.
ctest FunctionalTestJigsawApollo to validate this output. [#5710](https://github.com/DOI-USGS/ISIS3/issues/5710)
- Added OFFBODY and OFFBODYTRIM parameters to cam2cam. Added tests and updated documentation. [#3602] (https://github.com/DOI-USGS/ISIS3/issues/3602)
- Added support for reading, writing, and viewing GeoTIFFs in ISIS. [#5618](https://github.com/DOI-USGS/ISIS3/pull/5618)
- Added GDAL SRS propagation for systems outside of ISIS to display projected GTiffs. [#5736](https://github.com/DOI-USGS/ISIS3/pull/5736)

### Changed
- Enhanced csminit by removing the need to specify model and plugin [#5585](https://github.com/DOI-USGS/ISIS3/issues/5585)
Expand Down
2 changes: 1 addition & 1 deletion isis/src/base/objs/Blob/Blob.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -477,7 +477,7 @@ namespace Isis {
}

// update metadata
string jsonblobstr = pvl.toJson().dump();
string jsonblobstr = pvl.toJson()["Root"].dump();
string key = this->Type().toStdString() + "_" + this->Name().toStdString();
dataset->SetMetadataItem(key.c_str(), jsonblobstr.c_str(), "USGS");
}
Expand Down
176 changes: 106 additions & 70 deletions isis/src/base/objs/Cube/Cube.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1020,80 +1020,92 @@ namespace Isis {

isiscube.addObject(core);

if (dataset->GetSpatialRef()) {
char ** projStr = new char*[1];
const OGRSpatialReference &oSRS = *dataset->GetSpatialRef();
oSRS.exportToProj4(projStr);
QString qProjStr = QString::fromStdString(std::string(projStr[0]) + " +type=crs");
delete[] projStr[0];
delete[] projStr;

char ** projJsonStr = new char*[1];
oSRS.exportToPROJJSON(projJsonStr, nullptr);
nlohmann::json projJson = nlohmann::json::parse(projJsonStr[0]);
CPLFree(projJsonStr);

PvlGroup mappingGroup("Mapping");
mappingGroup.addKeyword(PvlKeyword("ProjectionName", "IProj"));
mappingGroup.addKeyword(PvlKeyword("EquatorialRadius", toString(oSRS.GetSemiMajor()), "meters"));
mappingGroup.addKeyword(PvlKeyword("PolarRadius", toString(oSRS.GetSemiMinor()), "meters"));

if (projJson.contains("base_crs")) {
projJson = projJson["base_crs"];
}
m_label->addObject(isiscube);
}

std::string direction = projJson["coordinate_system"]["axis"][1]["direction"];
if (direction == "east") {
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveEast"));
}
else if (direction == "west") {
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveWest"));
}
else {
QString msg = "Unknown direction [" + QString::fromStdString(direction) + "]";
throw IException(IException::Programmer, msg, _FILEINFO_);
}
mappingGroup.addKeyword(PvlKeyword("LongitudeDomain", "180"));
if (dataset->GetSpatialRef()) {
char ** projStr = new char*[1];
const OGRSpatialReference &oSRS = *dataset->GetSpatialRef();
oSRS.exportToProj4(projStr);
QString qProjStr = QString::fromStdString(std::string(projStr[0]) + " +type=crs");
delete[] projStr[0];
delete[] projStr;

char ** projJsonStr = new char*[1];
oSRS.exportToPROJJSON(projJsonStr, nullptr);
nlohmann::json projJson = nlohmann::json::parse(projJsonStr[0]);
CPLFree(projJsonStr);

PvlGroup mappingGroup("Mapping");
mappingGroup.addKeyword(PvlKeyword("ProjectionName", "IProj"));
mappingGroup.addKeyword(PvlKeyword("EquatorialRadius", toString(oSRS.GetSemiMajor()), "meters"));
mappingGroup.addKeyword(PvlKeyword("PolarRadius", toString(oSRS.GetSemiMinor()), "meters"));

if (projJson.contains("base_crs")) {
projJson = projJson["base_crs"];
}

std::string direction = projJson["coordinate_system"]["axis"][1]["direction"];
if (direction == "east") {
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveEast"));
}
else if (direction == "west") {
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveWest"));
}
else {
QString msg = "Unknown direction [" + QString::fromStdString(direction) + "]";
throw IException(IException::Programmer, msg, _FILEINFO_);
}

if (oSRS.GetSemiMajor() == oSRS.GetSemiMinor()) {
mappingGroup.addKeyword(PvlKeyword("LatitudeType", "Planetocentric"));
mappingGroup.addKeyword(PvlKeyword("ProjStr", qProjStr));

// Read the GeoTransform and get the elements we care about
double *padfTransform = new double[6];
dataset->GetGeoTransform(padfTransform);
if (abs(padfTransform[1]) != abs(padfTransform[5])) {
delete[] padfTransform;
QString msg = "Vertical and horizontal resolution do not match";
throw IException(IException::Io, msg, _FILEINFO_);
}
}
else {
mappingGroup.addKeyword(PvlKeyword("LatitudeType", "Planetographic"));
}

double dfScale;
double dfRes;
double upperLeftX;
double upperLeftY;
dfRes = padfTransform[1] * oSRS.GetLinearUnits();
upperLeftX = padfTransform[0];
upperLeftY = padfTransform[3];
if (oSRS.IsProjected()) {
const double dfDegToMeter = oSRS.GetSemiMajor() * M_PI / 180.0;
dfScale = dfDegToMeter / dfRes;
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "meters/pixel"));
}
else if (oSRS.IsGeographic()) {
dfScale = 1.0 / dfRes;
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "degrees/pixel"));
}
else {
QString msg = "Gdal spatial reference is not Geographic or Projected";
throw IException(IException::Io, msg, _FILEINFO_);
}
mappingGroup.addKeyword(PvlKeyword("Scale", toString(dfScale), "pixels/degree"));
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerX", toString(upperLeftX)));
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerY", toString(upperLeftY)));
mappingGroup.addKeyword(PvlKeyword("LongitudeDomain", "180"));
mappingGroup.addKeyword(PvlKeyword("ProjStr", qProjStr));

// Read the GeoTransform and get the elements we care about
double *padfTransform = new double[6];
dataset->GetGeoTransform(padfTransform);
if (abs(padfTransform[1]) != abs(padfTransform[5])) {
delete[] padfTransform;
QString msg = "Vertical and horizontal resolution do not match";
throw IException(IException::Io, msg, _FILEINFO_);
}

isiscube.addGroup(mappingGroup);
double dfScale;
double dfRes;
double upperLeftX;
double upperLeftY;
dfRes = padfTransform[1] * oSRS.GetLinearUnits();
upperLeftX = padfTransform[0];
upperLeftY = padfTransform[3];
if (oSRS.IsProjected()) {
const double dfDegToMeter = oSRS.GetSemiMajor() * M_PI / 180.0;
dfScale = dfDegToMeter / dfRes;
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "meters/pixel"));
}
m_label->addObject(isiscube);
else if (oSRS.IsGeographic()) {
dfScale = 1.0 / dfRes;
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "degrees/pixel"));
}
else {
QString msg = "Gdal spatial reference is not Geographic or Projected";
throw IException(IException::Io, msg, _FILEINFO_);
}
mappingGroup.addKeyword(PvlKeyword("Scale", toString(dfScale), "pixels/degree"));
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerX", toString(upperLeftX)));
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerY", toString(upperLeftY)));
delete[] padfTransform;

PvlObject &isiscube = m_label->findObject("IsisCube");
if (isiscube.hasGroup("Mapping")) {
isiscube.deleteGroup("Mapping");
}
isiscube.addGroup(mappingGroup);
}
GDALClose(dataset);

Expand Down Expand Up @@ -1122,7 +1134,7 @@ namespace Isis {
*
*/
void Cube::reopen(QString access) {
if (!m_labelFile) {
if (!isOpen()) {
QString msg = "Cube has not been opened yet. The filename to re-open is "
"unknown";
throw IException(IException::Programmer, msg, _FILEINFO_);
Expand Down Expand Up @@ -3024,7 +3036,6 @@ namespace Isis {
// update metadata
nlohmann::ordered_json jsonblob = this->label()->toJson()["Root"];
nlohmann::ordered_json jsonOut;
// std::cout << jsonblob << std::endl;
for (auto& [key, val] : jsonblob.items()) {
if (!val.contains("Bytes") || key == "Label") {
jsonOut[key] = val;
Expand All @@ -3033,6 +3044,31 @@ namespace Isis {
std::string jsonblobstr = jsonOut.dump();
std::string name = "CubeLabel";
gdalDataset()->SetMetadataItem(name.c_str(), jsonblobstr.c_str(), "USGS");

if (this->label()->findObject("IsisCube").hasGroup("Mapping")) {
PvlGroup &mappingGroup = this->label()->findObject("IsisCube").findGroup("Mapping");

if (mappingGroup.hasKeyword("ProjStr")) {
OGRSpatialReference *oSRS = new OGRSpatialReference();
oSRS->SetFromUserInput(mappingGroup.findKeyword("ProjStr")[0].toStdString().c_str());

gdalDataset()->SetSpatialRef(oSRS);

double dfRes = (double)mappingGroup.findKeyword("PixelResolution") / oSRS->GetLinearUnits();
double upperLeftX = (double) mappingGroup.findKeyword("UpperLeftCornerX");
double upperLeftY = (double) mappingGroup.findKeyword("UpperLeftCornerY");
double *padfTransform = new double[6];
padfTransform[1] = dfRes;
padfTransform[5] = -dfRes;
padfTransform[0] = upperLeftX;
padfTransform[3] = upperLeftY;
gdalDataset()->SetGeoTransform(padfTransform);
delete oSRS;
delete[] padfTransform;
}
}


return;
}

Expand Down
7 changes: 1 addition & 6 deletions isis/src/base/objs/Pvl/Pvl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,12 +155,7 @@ namespace Isis {
return jsonobj;
};

if (this->name() == "Root" && this->objects() == 1) {
return pvlobject_to_json(this->object(0));
}
else {
return pvlobject_to_json(*this);
}
return pvlobject_to_json(*this);
}

/**
Expand Down
74 changes: 71 additions & 3 deletions isis/tests/GTiffTests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ using json = nlohmann::json;
#include "TableRecord.h"

#include "CubeFixtures.h"
#include "TiffFixtures.h"
#include "TestUtilities.h"

#include "gmock/gmock.h"
Expand Down Expand Up @@ -174,9 +175,7 @@ TEST_P(GdalDnTypeGenerator, TestGTiffCreateWrite) {
in.close();
}

// Add Test for GTiff with no ISIS metadata

TEST_F(TempTestingFiles, TableTestsWriteReadGdal) {
TEST_F(TempTestingFiles, TestGTiffTableWriteRead) {
TableField f1("Column1", TableField::Integer);
TableField f2("Column2", TableField::Double);
TableField f3("Column3", TableField::Text, 10);
Expand Down Expand Up @@ -222,4 +221,73 @@ TEST_F(TempTestingFiles, TableTestsWriteReadGdal) {
for (int i = 0; i < t.Records(); i++) {
EXPECT_EQ(TableRecord::toString(t[i]).toStdString(), TableRecord::toString(t2[i]).toStdString());
}
}

TEST_F(ReadWriteTiff, TestGTiffSRS) {
PvlGroup mapping("Mapping");
mapping.addKeyword(PvlKeyword("ProjectionName", "IProj"));
mapping.addKeyword(PvlKeyword("EquatorialRadius", "3396190.0", "meters"));
mapping.addKeyword(PvlKeyword("PolarRadius", "3396190.0", "meters"));
mapping.addKeyword(PvlKeyword("LongitudeDirection", "PositiveEast"));
mapping.addKeyword(PvlKeyword("LatitudeType", "Planetocentric"));
mapping.addKeyword(PvlKeyword("LongitudeDomain", "180"));
mapping.addKeyword(PvlKeyword("ProjStr", "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs +type=crs"));
mapping.addKeyword(PvlKeyword("PixelResolution", "20.0", "meters/pixel"));
mapping.addKeyword(PvlKeyword("Scale", "2963.7348761653", "pixels/degree"));
mapping.addKeyword(PvlKeyword("UpperLeftCornerX", "8120050.0"));
mapping.addKeyword(PvlKeyword("UpperLeftCornerY", "-230230.0"));

createTiff(SignedWord);

Cube tiff;
tiff.open(path, "rw");
tiff.putGroup(mapping);
tiff.reopen();

Pvl &label = *tiff.label();
PvlGroup returnMapping = label.findObject("IsisCube").findGroup("Mapping");
EXPECT_TRUE(returnMapping == mapping);

dataset = GDALDataset::FromHandle(GDALOpen(path.toStdString().c_str(), GA_ReadOnly));

char ** projStr = new char*[1];
ASSERT_TRUE(dataset->GetSpatialRef() != nullptr);
const OGRSpatialReference &oSRS = *dataset->GetSpatialRef();
oSRS.exportToProj4(projStr);
std::string projAsCPPStr = projStr[0];
EXPECT_EQ(projAsCPPStr, "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs");

char ** projJsonStr = new char*[1];
oSRS.exportToPROJJSON(projJsonStr, nullptr);
nlohmann::ordered_json projJson = nlohmann::ordered_json::parse(projJsonStr[0]);
CPLFree(projJsonStr);

ASSERT_TRUE(projJson.contains("base_crs"));
ASSERT_TRUE(projJson["base_crs"].contains("coordinate_system"));
ASSERT_TRUE(projJson["base_crs"].contains("datum"));
ASSERT_TRUE(projJson["base_crs"]["coordinate_system"].contains("axis"));
ASSERT_TRUE(projJson["base_crs"]["datum"].contains("ellipsoid"));
ASSERT_TRUE(projJson["base_crs"]["datum"]["ellipsoid"].contains("radius"));

EXPECT_EQ(projJson["base_crs"]["coordinate_system"]["axis"][1]["direction"], "east");
EXPECT_EQ(projJson["base_crs"]["datum"]["ellipsoid"]["radius"], 3396190);

ASSERT_TRUE(projJson.contains("conversion"));
ASSERT_TRUE(projJson["conversion"].contains("name"));
ASSERT_TRUE(projJson["conversion"].contains("parameters"));
ASSERT_EQ(projJson["conversion"]["parameters"].size(), 5);

EXPECT_EQ(projJson["conversion"]["name"], "Equidistant Cylindrical");
EXPECT_EQ(projJson["conversion"]["parameters"][0]["name"], "Latitude of 1st standard parallel");
EXPECT_EQ(projJson["conversion"]["parameters"][0]["value"], 0);
EXPECT_EQ(projJson["conversion"]["parameters"][1]["name"], "Latitude of natural origin");
EXPECT_EQ(projJson["conversion"]["parameters"][1]["value"], 0);
EXPECT_EQ(projJson["conversion"]["parameters"][2]["name"], "Longitude of natural origin");
EXPECT_EQ(projJson["conversion"]["parameters"][2]["value"], 0);
EXPECT_EQ(projJson["conversion"]["parameters"][3]["name"], "False easting");
EXPECT_EQ(projJson["conversion"]["parameters"][3]["value"], 0);
EXPECT_EQ(projJson["conversion"]["parameters"][4]["name"], "False northing");
EXPECT_EQ(projJson["conversion"]["parameters"][4]["value"], 0);

dataset->Close();
}