Merge pull request #9 from NASA-IMPACT/fixDatum

add +datum=WGS84 to the crs

.github/workflows/build_tox_tests.yml

@@ -3,7 +3,7 @@ name: build_tox_tests
 on: [push]
 
 jobs:
-  build_push_to_ecr:
+  tox_in_docker:
     runs-on: ubuntu-latest
 
     steps:

hls_hdf_to_cog/hls_hdf_to_cog.py

@@ -32,7 +32,7 @@
 from rasterio.rio import options
 from rio_cogeo.cogeo import cog_translate, cog_validate
 from rio_cogeo.profiles import cog_profiles
-
+from rasterio.dtypes import _gdal_typename
 
 S30_BAND_NAMES = (
     "B01",
@@ -141,15 +141,43 @@ def main(input, output_dir, product, cogeo_profile, blocksize, creation_options)
                 output_name = os.path.join(output_dir, fname)
 
                 with rasterio.open(sds) as sub_dst:
-                    cog_translate(
-                        sub_dst,
-                        output_name,
-                        output_profile,
-                        config=config,
-                        forward_band_tags=True,
-                        overview_resampling="nearest",
-                        quiet=True,
-                    )
+                    # add the datum to the CRS
+                    # TODO: add a test
+                    new_crs = sub_dst.crs.to_proj4() + " +datum=WGS84"
+
+                    # We create a InMemory dataset using `GDAL MEM driver`
+                    # ref: https://github.com/rasterio/rasterio/blob/master/rasterio/_io.pyx#L1946-L1955
+                    data = sub_dst.read()
+                    info = {
+                        "DATAPOINTER": data.__array_interface__["data"][0],
+                        "PIXELS": sub_dst.width,
+                        "LINES": sub_dst.height,
+                        "BANDS": sub_dst.count,
+                        "DATATYPE": _gdal_typename(data.dtype.name),
+                    }
+                    dataset_options = ",".join(f"{name}={val}" for name, val in info.items())
+                    datasetname = f"MEM:::{dataset_options}"
+                    with rasterio.open(datasetname, "r+") as src:
+                        src.nodata = sub_dst.nodata  # set nodata
+                        src.crs = new_crs  # set CRS
+                        src.transform = sub_dst.transform  # add geotransform
+                        # set Metadata
+                        src.update_tags(**sub_dst.tags())
+                        src.colorinterp = sub_dst.colorinterp
+                        for i, b in enumerate(sub_dst.indexes):
+                            src.set_band_description(i + 1, sub_dst.descriptions[b - 1])
+                            src.update_tags(i + 1, **sub_dst.tags(b))
+
+                        # Create COG from the `Mem` dataset
+                        cog_translate(
+                            src,
+                            output_name,
+                            output_profile,
+                            config=config,
+                            forward_band_tags=True,
+                            overview_resampling="nearest",
+                            quiet=True,
+                        )
 
                 assert cog_validate(output_name)
 

tests/test_hls_hdf_to_cog.py

@@ -1,11 +1,13 @@
 import os
+import rasterio
 from click.testing import CliRunner
 from hls_hdf_to_cog.hls_hdf_to_cog import main
 
 
 current_dir = os.path.dirname(__file__)
 test_dir = os.path.join(current_dir, "data")
 data_dir = "/hls-testing_data"
+datum_string = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84"'
 
 
 def test_hls_hdf_to_cog_S30():
@@ -16,6 +18,9 @@ def test_hls_hdf_to_cog_S30():
                                   "--product", "S30"])
     print(result.exception)
     assert result.exit_code == 0
+    outputfile = os.path.join(test_dir, granule_basename.format(".B01.tif"))
+    with rasterio.open(outputfile, "r") as src:
+        assert datum_string in src.crs.to_wkt()
 
 
 def test_hls_hdf_to_cog_L30():
@@ -24,8 +29,11 @@ def test_hls_hdf_to_cog_L30():
     runner = CliRunner()
     result = runner.invoke(main, [inputfile, "--output-dir", test_dir,
                                   "--product", "L30"])
-    print(result.exception)
     assert result.exit_code == 0
+    print(result.exception)
+    outputfile = os.path.join(test_dir, granule_basename.format(".B01.tif"))
+    with rasterio.open(outputfile, "r") as src:
+        assert datum_string in src.crs.to_wkt()
 
 
 def test_hls_hdf_to_cog_S30_Angle():
@@ -36,3 +44,6 @@ def test_hls_hdf_to_cog_S30_Angle():
                                   "--product", "S30_ANGLES"])
     print(result.exception)
     assert result.exit_code == 0
+    outputfile = os.path.join(test_dir, "HLS.S30.T35JMG.2020192T074619.v1.5.SAA.tif")
+    with rasterio.open(outputfile, "r") as src:
+        assert datum_string in src.crs.to_wkt()