diff --git a/ale/base/base.py b/ale/base/base.py index bf6ad36fa..1df9c7e0b 100644 --- a/ale/base/base.py +++ b/ale/base/base.py @@ -1,6 +1,9 @@ import pvl import json +import tempfile +import os + class Driver(): """ Base class for all Drivers. @@ -323,3 +326,61 @@ def center_ephemeris_time(self): @property def short_mission_name(self): return self.__module__.split('.')[-1].split('_')[0] + + @property + def projection(self): + if not hasattr(self, "_projection"): + try: + from osgeo import gdal + except: + self._projection = "" + return self._projection + + if isinstance(self._file, pvl.PVLModule): + # save it to a temp folder + with tempfile.NamedTemporaryFile() as tmp: + tmp.write(pvl.dumps(self._file)) + + geodata = gdal.Open(tempfile.name) + self._projection = geodata.GetSpatialRef() + else: + # should be a path + if not os.path.exists(self._file): + self._projection = "" + else: + geodata = gdal.Open(self._file) + self._projection = geodata.GetSpatialRef() + + # is None if not projected + if self._projection: + self._projection = self._projection.ExportToProj4() + else: + self._projection = "" + return self._projection + + + @property + def geotransform(self): + if not hasattr(self, "_geotransform"): + try: + from osgeo import gdal + except: + self._geotransform = (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) + return self._geotransform + + if isinstance(self._file, pvl.PVLModule): + # save it to a temp folder + with tempfile.NamedTemporaryFile() as tmp: + tmp.write(pvl.dumps(self._file)) + + geodata = gdal.Open(tempfile.name) + self._geotransform = geodata.GetGeoTransform() + else: + # should be a path + if not os.path.exists(self._file): + self._geotransform = (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) + else: + geodata = gdal.Open(self._file) + self._geotransform = geodata.GetGeoTransform() + + return self._geotransform \ No newline at end of file diff --git a/ale/formatters/formatter.py b/ale/formatters/formatter.py index b33287bbf..23968d627 100644 --- a/ale/formatters/formatter.py +++ b/ale/formatters/formatter.py @@ -196,6 +196,8 @@ def to_isd(driver): meta_data['sun_position'] = sun_position + meta_data["projection"] = driver.projection + meta_data["geotransform"] = driver.geotransform # check that there is a valid sensor model name if 'name_model' not in meta_data: diff --git a/ale/formatters/usgscsm_formatter.py b/ale/formatters/usgscsm_formatter.py index d7fd7dba7..44438fdac 100644 --- a/ale/formatters/usgscsm_formatter.py +++ b/ale/formatters/usgscsm_formatter.py @@ -51,6 +51,9 @@ def to_usgscsm(driver): 'unit' : 'm' } + isd_data["projection"] = driver.projection + isd_data["geotransform"] = driver.geotransform + # shared isd keywords for Framer and Linescanner if isinstance(driver, LineScanner) or isinstance(driver, Framer): # exterior orientation for just Framer and LineScanner diff --git a/environment.yml b/environment.yml index 7b84927f7..b4c11109f 100644 --- a/environment.yml +++ b/environment.yml @@ -7,6 +7,7 @@ dependencies: - cmake>=3.15 - pytest - eigen + - gdal - jupyter - nlohmann_json - numpy diff --git a/include/ale/Util.h b/include/ale/Util.h index d27244f30..dbc50be69 100644 --- a/include/ale/Util.h +++ b/include/ale/Util.h @@ -34,6 +34,8 @@ namespace ale { std::string getPlatformName(nlohmann::json isd); std::string getLogFile(nlohmann::json isd); std::string getIsisCameraVersion(nlohmann::json isd); + std::string getProjection(nlohmann::json isd); + int getTotalLines(nlohmann::json isd); int getTotalSamples(nlohmann::json isd); double getStartingTime(nlohmann::json isd); @@ -45,6 +47,7 @@ namespace ale { double getFocalLengthUncertainty(nlohmann::json isd); std::vector getFocal2PixelLines(nlohmann::json isd); std::vector getFocal2PixelSamples(nlohmann::json isd); + std::vector getGeoTransform(nlohmann::json isd); double getDetectorCenterLine(nlohmann::json isd); double getDetectorCenterSample(nlohmann::json isd); double getDetectorStartingLine(nlohmann::json isd); diff --git a/recipe/meta.yaml b/recipe/meta.yaml index c96b750a3..7ca81bb52 100644 --- a/recipe/meta.yaml +++ b/recipe/meta.yaml @@ -31,7 +31,9 @@ requirements: - scipy >=1.4.0 - spiceypy >=4.0.1 - pyyaml - + run_contrained: + - gdal + test: imports: - ale diff --git a/src/Util.cpp b/src/Util.cpp index 297c46c1c..25fbd401f 100644 --- a/src/Util.cpp +++ b/src/Util.cpp @@ -37,6 +37,31 @@ std::string getImageId(json isd) { return id; } + +std::vector getGeoTransform(json isd) { + std::vector transform = {}; + try { + transform = isd.at("geotransform").get>(); + } catch (std::exception &e) { + std::string originalError = e.what(); + std::string msg = "Could not parse the geo_transform. ERROR: \n" + originalError; + throw std::runtime_error(msg); + } + return transform; +} + + +std::string getProjection(json isd) { + std::string projection_string = ""; + try { + projection_string = isd.at("projection"); + } catch (...) { + throw std::runtime_error("Could not parse the projection string."); + } + return projection_string; +} + + std::string getSensorName(json isd) { std::string name = ""; try { diff --git a/tests/ctests/IsdTests.cpp b/tests/ctests/IsdTests.cpp index 23fe47737..fb6a34dc6 100644 --- a/tests/ctests/IsdTests.cpp +++ b/tests/ctests/IsdTests.cpp @@ -77,6 +77,22 @@ TEST(Isd, LogFile) { EXPECT_STREQ(ale::getLogFile(j).c_str(), "fake/path"); } + +TEST(Isd, Projection) { + nlohmann::json proj; + proj["projection"] = "totally a proj4 string"; + proj["geotransform"] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; + + std::vector coeffs = ale::getGeoTransform(proj); + std::string proj4str = ale::getProjection(proj); + EXPECT_EQ(proj4str, "totally a proj4 string"); + ASSERT_EQ(coeffs.size(), 6); + EXPECT_DOUBLE_EQ(coeffs[1], 1.0); + EXPECT_DOUBLE_EQ(coeffs[5], 1.0); + EXPECT_DOUBLE_EQ(coeffs[0], 0.0); +} + + TEST(Isd, TransverseDistortion) { nlohmann::json trans; trans["optical_distortion"]["transverse"]["x"] = {1}; @@ -571,6 +587,21 @@ TEST(Isd, BadStartingDetector) { } } + +TEST(Isd, BadProjection) { + nlohmann::json proj; + proj["projection"] = NULL; + + try { + std::string proj4str = ale::getProjection(proj); + FAIL() << "Expected exception to be thrown"; + } + catch(std::exception &e) { + EXPECT_EQ(std::string(e.what()), "Could not parse the projection string."); + } +} + + TEST(Isd, BadFocal2Pixel) { std::string bad_json_str("{}"); try { diff --git a/tests/pytests/data/B10_013341_1010_XN_79S172W/B10_013341_1010_XN_79S172W_isis3.lbl b/tests/pytests/data/B10_013341_1010_XN_79S172W/B10_013341_1010_XN_79S172W_isis3.lbl index 111547542..504171f00 100644 --- a/tests/pytests/data/B10_013341_1010_XN_79S172W/B10_013341_1010_XN_79S172W_isis3.lbl +++ b/tests/pytests/data/B10_013341_1010_XN_79S172W/B10_013341_1010_XN_79S172W_isis3.lbl @@ -69,6 +69,25 @@ Object = IsisCube Source = isis End_Group + Group = Mapping + ProjectionName = Sinusoidal + CenterLongitude = 148.36859083039 + TargetName = MARS + EquatorialRadius = 3396190.0 + PolarRadius = 3376200.0 + LatitudeType = Planetocentric + LongitudeDirection = PositiveEast + LongitudeDomain = 360 + MinimumLatitude = 63.636322793577 + MaximumLatitude = 87.296295823424 + MinimumLongitude = 139.6658284858 + MaximumLongitude = 157.07135317498 + UpperLeftCornerX = -219771.1526456 + UpperLeftCornerY = 5175537.8728989 + PixelResolution = 1455.4380969907 + Scale = 40.726361118253 + End_Group + Group = AlphaCube AlphaSamples = 5000 AlphaLines = 24576 diff --git a/tests/pytests/test_load.py b/tests/pytests/test_load.py index 20776a9f2..8c7e291c3 100644 --- a/tests/pytests/test_load.py +++ b/tests/pytests/test_load.py @@ -73,12 +73,10 @@ def test_load_mes_from_metakernels(tmpdir, monkeypatch, mess_kernels): mk_file.write(mk_str) usgscsm_isd_obj = ale.load(label_file, verbose=True) - assert usgscsm_isd_obj['name_platform'] == 'MESSENGER' assert usgscsm_isd_obj['name_sensor'] == 'MERCURY DUAL IMAGING SYSTEM NARROW ANGLE CAMERA' assert usgscsm_isd_obj['name_model'] == 'USGS_ASTRO_FRAME_SENSOR_MODEL' - def test_load_mes_with_no_metakernels(tmpdir, monkeypatch, mess_kernels): monkeypatch.setenv('ALESPICEROOT', str(tmpdir)) diff --git a/tests/pytests/test_mex_drivers.py b/tests/pytests/test_mex_drivers.py index 650d9b478..ec5035c28 100644 --- a/tests/pytests/test_mex_drivers.py +++ b/tests/pytests/test_mex_drivers.py @@ -202,7 +202,9 @@ def usgscsm_compare_dict(): "t0_ephemeris": -101.83713859319687, "dt_ephemeris": 40.734855437278746, "t0_quaternion": -101.83713859319687, - "dt_quaternion": 40.734855437278746 + "dt_quaternion": 40.734855437278746, + "projection" : "", + "geotransform" : (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) }, "isis" : diff --git a/tests/pytests/test_usgscsm_formatter.py b/tests/pytests/test_usgscsm_formatter.py index 9a9ae6d60..70ecb641c 100644 --- a/tests/pytests/test_usgscsm_formatter.py +++ b/tests/pytests/test_usgscsm_formatter.py @@ -9,6 +9,8 @@ from ale.base.data_naif import NaifSpice from ale.rotation import ConstantRotation, TimeDependentRotation +from conftest import get_image_label + class TestDriver(Driver, NaifSpice): """ Test Driver implementation with dummy values @@ -342,3 +344,20 @@ def test_line_scan_sun_position(test_line_scan_driver): assert sun_position_obj['positions'] == [[0, 1, 2], [3, 4, 5]] assert sun_position_obj['velocities'] == [[0, -1, -2], [-3, -4, -5]] assert sun_position_obj['unit'] == 'm' + +def test_no_projection(test_frame_driver): + isd = usgscsm_formatter.to_usgscsm(test_frame_driver) + # isn't using real projection so it should be None + assert isd['projection'] == None + +def test_isis_projection(): + isd = usgscsm_formatter.to_usgscsm(TestLineScanner(get_image_label('B10_013341_1010_XN_79S172W', "isis3"))) + assert isd["projection"] == "+proj=sinu +lon_0=148.36859083039 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs" + + +def test_isis_geotransform(): + isd = usgscsm_formatter.to_usgscsm(TestLineScanner(get_image_label('B10_013341_1010_XN_79S172W', "isis3"))) + expected = (-219771.1526456, 1455.4380969907, 0.0, 5175537.8728989, 0.0, -1455.4380969907) + for value, truth in zip(isd["geotransform"], expected): + pytest.approx(value, truth) +