Skip to content

Commit 463db72

Browse files
acpaquetteacpaquette
authored
GTiff SRS Propagation (#5736)
* Set the SpatialReference for GTiffs that use the proj projection plugin * Added changelog entry * Added test for opening a non-gdal created ISIS cube and writing a mapping group to it * Fixed opening/reopening with blobs * Fixed GTiff test name
1 parent 1e28c61 commit 463db72

File tree

5 files changed

+180
-80
lines changed

5 files changed

+180
-80
lines changed

‎CHANGELOG.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ release.
4040
ctest FunctionalTestJigsawApollo to validate this output. [#5710](https://github.com/DOI-USGS/ISIS3/issues/5710)
4141
- Added OFFBODY and OFFBODYTRIM parameters to cam2cam. Added tests and updated documentation. [#3602] (https://github.com/DOI-USGS/ISIS3/issues/3602)
4242
- Added support for reading, writing, and viewing GeoTIFFs in ISIS. [#5618](https://github.com/DOI-USGS/ISIS3/pull/5618)
43+
- Added GDAL SRS propagation for systems outside of ISIS to display projected GTiffs. [#5736](https://github.com/DOI-USGS/ISIS3/pull/5736)
4344

4445
### Changed
4546
- Enhanced csminit by removing the need to specify model and plugin [#5585](https://github.com/DOI-USGS/ISIS3/issues/5585)

‎isis/src/base/objs/Blob/Blob.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -477,7 +477,7 @@ namespace Isis {
477477
}
478478

479479
// update metadata
480-
string jsonblobstr = pvl.toJson().dump();
480+
string jsonblobstr = pvl.toJson()["Root"].dump();
481481
string key = this->Type().toStdString() + "_" + this->Name().toStdString();
482482
dataset->SetMetadataItem(key.c_str(), jsonblobstr.c_str(), "USGS");
483483
}

‎isis/src/base/objs/Cube/Cube.cpp

Lines changed: 106 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -1020,80 +1020,92 @@ namespace Isis {
10201020

10211021
isiscube.addObject(core);
10221022

1023-
if (dataset->GetSpatialRef()) {
1024-
char ** projStr = new char*[1];
1025-
const OGRSpatialReference &oSRS = *dataset->GetSpatialRef();
1026-
oSRS.exportToProj4(projStr);
1027-
QString qProjStr = QString::fromStdString(std::string(projStr[0]) + " +type=crs");
1028-
delete[] projStr[0];
1029-
delete[] projStr;
1030-
1031-
char ** projJsonStr = new char*[1];
1032-
oSRS.exportToPROJJSON(projJsonStr, nullptr);
1033-
nlohmann::json projJson = nlohmann::json::parse(projJsonStr[0]);
1034-
CPLFree(projJsonStr);
1035-
1036-
PvlGroup mappingGroup("Mapping");
1037-
mappingGroup.addKeyword(PvlKeyword("ProjectionName", "IProj"));
1038-
mappingGroup.addKeyword(PvlKeyword("EquatorialRadius", toString(oSRS.GetSemiMajor()), "meters"));
1039-
mappingGroup.addKeyword(PvlKeyword("PolarRadius", toString(oSRS.GetSemiMinor()), "meters"));
1040-
1041-
if (projJson.contains("base_crs")) {
1042-
projJson = projJson["base_crs"];
1043-
}
1023+
m_label->addObject(isiscube);
1024+
}
10441025

1045-
std::string direction = projJson["coordinate_system"]["axis"][1]["direction"];
1046-
if (direction == "east") {
1047-
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveEast"));
1048-
}
1049-
else if (direction == "west") {
1050-
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveWest"));
1051-
}
1052-
else {
1053-
QString msg = "Unknown direction [" + QString::fromStdString(direction) + "]";
1054-
throw IException(IException::Programmer, msg, _FILEINFO_);
1055-
}
1056-
mappingGroup.addKeyword(PvlKeyword("LongitudeDomain", "180"));
1026+
if (dataset->GetSpatialRef()) {
1027+
char ** projStr = new char*[1];
1028+
const OGRSpatialReference &oSRS = *dataset->GetSpatialRef();
1029+
oSRS.exportToProj4(projStr);
1030+
QString qProjStr = QString::fromStdString(std::string(projStr[0]) + " +type=crs");
1031+
delete[] projStr[0];
1032+
delete[] projStr;
1033+
1034+
char ** projJsonStr = new char*[1];
1035+
oSRS.exportToPROJJSON(projJsonStr, nullptr);
1036+
nlohmann::json projJson = nlohmann::json::parse(projJsonStr[0]);
1037+
CPLFree(projJsonStr);
1038+
1039+
PvlGroup mappingGroup("Mapping");
1040+
mappingGroup.addKeyword(PvlKeyword("ProjectionName", "IProj"));
1041+
mappingGroup.addKeyword(PvlKeyword("EquatorialRadius", toString(oSRS.GetSemiMajor()), "meters"));
1042+
mappingGroup.addKeyword(PvlKeyword("PolarRadius", toString(oSRS.GetSemiMinor()), "meters"));
1043+
1044+
if (projJson.contains("base_crs")) {
1045+
projJson = projJson["base_crs"];
1046+
}
1047+
1048+
std::string direction = projJson["coordinate_system"]["axis"][1]["direction"];
1049+
if (direction == "east") {
1050+
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveEast"));
1051+
}
1052+
else if (direction == "west") {
1053+
mappingGroup.addKeyword(PvlKeyword("LongitudeDirection", "PositiveWest"));
1054+
}
1055+
else {
1056+
QString msg = "Unknown direction [" + QString::fromStdString(direction) + "]";
1057+
throw IException(IException::Programmer, msg, _FILEINFO_);
1058+
}
1059+
1060+
if (oSRS.GetSemiMajor() == oSRS.GetSemiMinor()) {
10571061
mappingGroup.addKeyword(PvlKeyword("LatitudeType", "Planetocentric"));
1058-
mappingGroup.addKeyword(PvlKeyword("ProjStr", qProjStr));
1059-
1060-
// Read the GeoTransform and get the elements we care about
1061-
double *padfTransform = new double[6];
1062-
dataset->GetGeoTransform(padfTransform);
1063-
if (abs(padfTransform[1]) != abs(padfTransform[5])) {
1064-
delete[] padfTransform;
1065-
QString msg = "Vertical and horizontal resolution do not match";
1066-
throw IException(IException::Io, msg, _FILEINFO_);
1067-
}
1062+
}
1063+
else {
1064+
mappingGroup.addKeyword(PvlKeyword("LatitudeType", "Planetographic"));
1065+
}
10681066

1069-
double dfScale;
1070-
double dfRes;
1071-
double upperLeftX;
1072-
double upperLeftY;
1073-
dfRes = padfTransform[1] * oSRS.GetLinearUnits();
1074-
upperLeftX = padfTransform[0];
1075-
upperLeftY = padfTransform[3];
1076-
if (oSRS.IsProjected()) {
1077-
const double dfDegToMeter = oSRS.GetSemiMajor() * M_PI / 180.0;
1078-
dfScale = dfDegToMeter / dfRes;
1079-
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "meters/pixel"));
1080-
}
1081-
else if (oSRS.IsGeographic()) {
1082-
dfScale = 1.0 / dfRes;
1083-
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "degrees/pixel"));
1084-
}
1085-
else {
1086-
QString msg = "Gdal spatial reference is not Geographic or Projected";
1087-
throw IException(IException::Io, msg, _FILEINFO_);
1088-
}
1089-
mappingGroup.addKeyword(PvlKeyword("Scale", toString(dfScale), "pixels/degree"));
1090-
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerX", toString(upperLeftX)));
1091-
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerY", toString(upperLeftY)));
1067+
mappingGroup.addKeyword(PvlKeyword("LongitudeDomain", "180"));
1068+
mappingGroup.addKeyword(PvlKeyword("ProjStr", qProjStr));
1069+
1070+
// Read the GeoTransform and get the elements we care about
1071+
double *padfTransform = new double[6];
1072+
dataset->GetGeoTransform(padfTransform);
1073+
if (abs(padfTransform[1]) != abs(padfTransform[5])) {
10921074
delete[] padfTransform;
1075+
QString msg = "Vertical and horizontal resolution do not match";
1076+
throw IException(IException::Io, msg, _FILEINFO_);
1077+
}
10931078

1094-
isiscube.addGroup(mappingGroup);
1079+
double dfScale;
1080+
double dfRes;
1081+
double upperLeftX;
1082+
double upperLeftY;
1083+
dfRes = padfTransform[1] * oSRS.GetLinearUnits();
1084+
upperLeftX = padfTransform[0];
1085+
upperLeftY = padfTransform[3];
1086+
if (oSRS.IsProjected()) {
1087+
const double dfDegToMeter = oSRS.GetSemiMajor() * M_PI / 180.0;
1088+
dfScale = dfDegToMeter / dfRes;
1089+
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "meters/pixel"));
10951090
}
1096-
m_label->addObject(isiscube);
1091+
else if (oSRS.IsGeographic()) {
1092+
dfScale = 1.0 / dfRes;
1093+
mappingGroup.addKeyword(PvlKeyword("PixelResolution", toString(dfRes), "degrees/pixel"));
1094+
}
1095+
else {
1096+
QString msg = "Gdal spatial reference is not Geographic or Projected";
1097+
throw IException(IException::Io, msg, _FILEINFO_);
1098+
}
1099+
mappingGroup.addKeyword(PvlKeyword("Scale", toString(dfScale), "pixels/degree"));
1100+
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerX", toString(upperLeftX)));
1101+
mappingGroup.addKeyword(PvlKeyword("UpperLeftCornerY", toString(upperLeftY)));
1102+
delete[] padfTransform;
1103+
1104+
PvlObject &isiscube = m_label->findObject("IsisCube");
1105+
if (isiscube.hasGroup("Mapping")) {
1106+
isiscube.deleteGroup("Mapping");
1107+
}
1108+
isiscube.addGroup(mappingGroup);
10971109
}
10981110
GDALClose(dataset);
10991111

@@ -1122,7 +1134,7 @@ namespace Isis {
11221134
*
11231135
*/
11241136
void Cube::reopen(QString access) {
1125-
if (!m_labelFile) {
1137+
if (!isOpen()) {
11261138
QString msg = "Cube has not been opened yet. The filename to re-open is "
11271139
"unknown";
11281140
throw IException(IException::Programmer, msg, _FILEINFO_);
@@ -3024,7 +3036,6 @@ namespace Isis {
30243036
// update metadata
30253037
nlohmann::ordered_json jsonblob = this->label()->toJson()["Root"];
30263038
nlohmann::ordered_json jsonOut;
3027-
// std::cout << jsonblob << std::endl;
30283039
for (auto& [key, val] : jsonblob.items()) {
30293040
if (!val.contains("Bytes") || key == "Label") {
30303041
jsonOut[key] = val;
@@ -3033,6 +3044,31 @@ namespace Isis {
30333044
std::string jsonblobstr = jsonOut.dump();
30343045
std::string name = "CubeLabel";
30353046
gdalDataset()->SetMetadataItem(name.c_str(), jsonblobstr.c_str(), "USGS");
3047+
3048+
if (this->label()->findObject("IsisCube").hasGroup("Mapping")) {
3049+
PvlGroup &mappingGroup = this->label()->findObject("IsisCube").findGroup("Mapping");
3050+
3051+
if (mappingGroup.hasKeyword("ProjStr")) {
3052+
OGRSpatialReference *oSRS = new OGRSpatialReference();
3053+
oSRS->SetFromUserInput(mappingGroup.findKeyword("ProjStr")[0].toStdString().c_str());
3054+
3055+
gdalDataset()->SetSpatialRef(oSRS);
3056+
3057+
double dfRes = (double)mappingGroup.findKeyword("PixelResolution") / oSRS->GetLinearUnits();
3058+
double upperLeftX = (double) mappingGroup.findKeyword("UpperLeftCornerX");
3059+
double upperLeftY = (double) mappingGroup.findKeyword("UpperLeftCornerY");
3060+
double *padfTransform = new double[6];
3061+
padfTransform[1] = dfRes;
3062+
padfTransform[5] = -dfRes;
3063+
padfTransform[0] = upperLeftX;
3064+
padfTransform[3] = upperLeftY;
3065+
gdalDataset()->SetGeoTransform(padfTransform);
3066+
delete oSRS;
3067+
delete[] padfTransform;
3068+
}
3069+
}
3070+
3071+
30363072
return;
30373073
}
30383074

‎isis/src/base/objs/Pvl/Pvl.cpp

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -155,12 +155,7 @@ namespace Isis {
155155
return jsonobj;
156156
};
157157

158-
if (this->name() == "Root" && this->objects() == 1) {
159-
return pvlobject_to_json(this->object(0));
160-
}
161-
else {
162-
return pvlobject_to_json(*this);
163-
}
158+
return pvlobject_to_json(*this);
164159
}
165160

166161
/**

‎isis/tests/GTiffTests.cpp

Lines changed: 71 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ using json = nlohmann::json;
1818
#include "TableRecord.h"
1919

2020
#include "CubeFixtures.h"
21+
#include "TiffFixtures.h"
2122
#include "TestUtilities.h"
2223

2324
#include "gmock/gmock.h"
@@ -174,9 +175,7 @@ TEST_P(GdalDnTypeGenerator, TestGTiffCreateWrite) {
174175
in.close();
175176
}
176177

177-
// Add Test for GTiff with no ISIS metadata
178-
179-
TEST_F(TempTestingFiles, TableTestsWriteReadGdal) {
178+
TEST_F(TempTestingFiles, TestGTiffTableWriteRead) {
180179
TableField f1("Column1", TableField::Integer);
181180
TableField f2("Column2", TableField::Double);
182181
TableField f3("Column3", TableField::Text, 10);
@@ -222,4 +221,73 @@ TEST_F(TempTestingFiles, TableTestsWriteReadGdal) {
222221
for (int i = 0; i < t.Records(); i++) {
223222
EXPECT_EQ(TableRecord::toString(t[i]).toStdString(), TableRecord::toString(t2[i]).toStdString());
224223
}
224+
}
225+
226+
TEST_F(ReadWriteTiff, TestGTiffSRS) {
227+
PvlGroup mapping("Mapping");
228+
mapping.addKeyword(PvlKeyword("ProjectionName", "IProj"));
229+
mapping.addKeyword(PvlKeyword("EquatorialRadius", "3396190.0", "meters"));
230+
mapping.addKeyword(PvlKeyword("PolarRadius", "3396190.0", "meters"));
231+
mapping.addKeyword(PvlKeyword("LongitudeDirection", "PositiveEast"));
232+
mapping.addKeyword(PvlKeyword("LatitudeType", "Planetocentric"));
233+
mapping.addKeyword(PvlKeyword("LongitudeDomain", "180"));
234+
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"));
235+
mapping.addKeyword(PvlKeyword("PixelResolution", "20.0", "meters/pixel"));
236+
mapping.addKeyword(PvlKeyword("Scale", "2963.7348761653", "pixels/degree"));
237+
mapping.addKeyword(PvlKeyword("UpperLeftCornerX", "8120050.0"));
238+
mapping.addKeyword(PvlKeyword("UpperLeftCornerY", "-230230.0"));
239+
240+
createTiff(SignedWord);
241+
242+
Cube tiff;
243+
tiff.open(path, "rw");
244+
tiff.putGroup(mapping);
245+
tiff.reopen();
246+
247+
Pvl &label = *tiff.label();
248+
PvlGroup returnMapping = label.findObject("IsisCube").findGroup("Mapping");
249+
EXPECT_TRUE(returnMapping == mapping);
250+
251+
dataset = GDALDataset::FromHandle(GDALOpen(path.toStdString().c_str(), GA_ReadOnly));
252+
253+
char ** projStr = new char*[1];
254+
ASSERT_TRUE(dataset->GetSpatialRef() != nullptr);
255+
const OGRSpatialReference &oSRS = *dataset->GetSpatialRef();
256+
oSRS.exportToProj4(projStr);
257+
std::string projAsCPPStr = projStr[0];
258+
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");
259+
260+
char ** projJsonStr = new char*[1];
261+
oSRS.exportToPROJJSON(projJsonStr, nullptr);
262+
nlohmann::ordered_json projJson = nlohmann::ordered_json::parse(projJsonStr[0]);
263+
CPLFree(projJsonStr);
264+
265+
ASSERT_TRUE(projJson.contains("base_crs"));
266+
ASSERT_TRUE(projJson["base_crs"].contains("coordinate_system"));
267+
ASSERT_TRUE(projJson["base_crs"].contains("datum"));
268+
ASSERT_TRUE(projJson["base_crs"]["coordinate_system"].contains("axis"));
269+
ASSERT_TRUE(projJson["base_crs"]["datum"].contains("ellipsoid"));
270+
ASSERT_TRUE(projJson["base_crs"]["datum"]["ellipsoid"].contains("radius"));
271+
272+
EXPECT_EQ(projJson["base_crs"]["coordinate_system"]["axis"][1]["direction"], "east");
273+
EXPECT_EQ(projJson["base_crs"]["datum"]["ellipsoid"]["radius"], 3396190);
274+
275+
ASSERT_TRUE(projJson.contains("conversion"));
276+
ASSERT_TRUE(projJson["conversion"].contains("name"));
277+
ASSERT_TRUE(projJson["conversion"].contains("parameters"));
278+
ASSERT_EQ(projJson["conversion"]["parameters"].size(), 5);
279+
280+
EXPECT_EQ(projJson["conversion"]["name"], "Equidistant Cylindrical");
281+
EXPECT_EQ(projJson["conversion"]["parameters"][0]["name"], "Latitude of 1st standard parallel");
282+
EXPECT_EQ(projJson["conversion"]["parameters"][0]["value"], 0);
283+
EXPECT_EQ(projJson["conversion"]["parameters"][1]["name"], "Latitude of natural origin");
284+
EXPECT_EQ(projJson["conversion"]["parameters"][1]["value"], 0);
285+
EXPECT_EQ(projJson["conversion"]["parameters"][2]["name"], "Longitude of natural origin");
286+
EXPECT_EQ(projJson["conversion"]["parameters"][2]["value"], 0);
287+
EXPECT_EQ(projJson["conversion"]["parameters"][3]["name"], "False easting");
288+
EXPECT_EQ(projJson["conversion"]["parameters"][3]["value"], 0);
289+
EXPECT_EQ(projJson["conversion"]["parameters"][4]["name"], "False northing");
290+
EXPECT_EQ(projJson["conversion"]["parameters"][4]["value"], 0);
291+
292+
dataset->Close();
225293
}

0 commit comments

Comments
 (0)