correclty compute burst byte offests

README.md

@@ -67,19 +67,18 @@ A full accounting of omitted datasets and differing fields can be found below:
 
 * Annotation
     * Noise
-        * No intentional omissions or nulls.
+        * No intentional omissions or deviations.
     * Calibration
-        * No intentional omissions or nulls.
+        * No intentional omissions or deviations
     * RFI
-        * No intentional omissions or nulls.
+        * No intentional omissions or deviations.
     * Product
         * `generalAnnotation/productInformation/platformHeading`
             - calculated as average of input Level 1 SLCs, not recalculated from Level-0 slices.
         * `imageAnnotation/imageInformation/azimuthPixelSpacing`
             - calculated as average of input Level 1 SLCs, not Level-0 slices.
         * `imageAnnotation/imageInformation/imageStatistics/outputDataMean` / `outputDataStdDev`
             - calculated using `np.mean`/`np.std` on valid data.
-        * `swathTiming/burstList/burst/byteOffset` components are set to `''`.
 * Measurement GeoTIFFs
     * Invalid data as denoted by `swathTiming/burstList/burst/firstValidSample` and `lastValidSample` are set to zero. This done by the ASF extractor, not this tool.
     * TIFF tags **that are not GeoTIFF tags** are omitted. See Product Specification Table 3-8 for full list.

src/burst2safe/measurement.py

@@ -5,6 +5,7 @@
 
 import numpy as np
 from osgeo import gdal, osr
+from tifffile import TiffFile
 
 from burst2safe.base import create_content_unit, create_data_object
 from burst2safe.product import GeoPoint
@@ -36,6 +37,7 @@ def __init__(self, burst_infos: Iterable[BurstInfo], gcps: Iterable[GeoPoint], i
         self.data_std = None
         self.size_bytes = None
         self.md5 = None
+        self.byte_offsets = []
 
     def get_data(self, band: int = 1) -> np.ndarray:
         """Get the data from the measurement from ASF burst GeoTIFFs.
@@ -53,6 +55,23 @@ def get_data(self, band: int = 1) -> np.ndarray:
             ds = None
         return data
 
+    def get_burst_byte_offsets(self):
+        with TiffFile(self.path) as tif:
+            if len(tif.pages) != 1:
+                raise ValueError('Byte offset calculation only valid for GeoTIFFs with one band.')
+            page = tif.pages[0]
+
+            if str(page.compression) != 'COMPRESSION.NONE':
+                raise ValueError('Byte offset calculation only valid for uncompressed GeoTIFFs.')
+
+            if page.chunks != (1, page.imagewidth):
+                raise ValueError('Byte offset calculation only valid for GeoTIFFs with one line per block.')
+
+            offsets = page.dataoffsets
+
+        byte_offsets = [offsets[self.burst_length * i] for i in range(len(self.burst_infos))]
+        return byte_offsets
+
     def add_metadata(self, dataset: gdal.Dataset):
         """Add metadata to an existing GDAL dataset.
 
@@ -86,10 +105,13 @@ def create_geotiff(self, out_path: Path, update_info=True):
         self.add_metadata(mem_ds)
         gdal.Translate(str(out_path), mem_ds, format='GTiff')
         mem_ds = None
-        self.path = out_path
+
         if update_info:
+            self.path = out_path
             self.data_mean = np.mean(data[data != 0])
             self.data_std = np.std(data[data != 0])
+            self.byte_offsets = self.get_burst_byte_offsets()
+
             with open(self.path, 'rb') as f:
                 file_bytes = f.read()
                 self.size_bytes = len(file_bytes)

src/burst2safe/product.py

@@ -215,6 +215,20 @@ def update_gcps(self):
             )
             self.gcps.append(gcp)
 
+    def update_burst_byte_offsets(self, byte_offsets: Iterable[int]):
+        """Update the byte offsets in the burstList element.
+
+        Args:
+            byte_offsets: The byte offsets to update
+        """
+        if self.swath_timing is None or self.xml.find('swathTiming') is None:
+            raise ValueError('Product must be assembled before updating burst byte offsets.')
+
+        for swath_timing in [self.swath_timing, self.xml.find('swathTiming')]:
+            burst_list = swath_timing.find('burstList')
+            for i, byte_offset in enumerate(byte_offsets):
+                burst_list[i].find('byteOffset').text = str(byte_offset)
+
     def create_geolocation_grid(self):
         """Create the geolocationGrid element."""
         geolocation_grid = ET.Element('geolocationGrid')

src/burst2safe/swath.py

@@ -129,6 +129,7 @@ def write(self, update_info: bool = True):
         """
         self.measurement.write(self.measurement_name)
         self.product.update_data_stats(self.measurement.data_mean, self.measurement.data_std)
+        self.product.update_burst_byte_offsets(self.measurement.byte_offsets)
         self.product.write(self.product_name)
         self.noise.write(self.noise_name)
         self.calibration.write(self.calibration_name)

tests/test_measurement.py

@@ -19,8 +19,8 @@ def burst_datas(burst_infos, tmp_path):
         n = i + 1
         burst_data = deepcopy(burst_info)
         burst_data.data_path = tmp_path / f'burst{n}.tif'
-        burst_data.length = 10
-        burst_data.width = 20
+        burst_data.length = 1000
+        burst_data.width = 2000
 
         shape = (burst_data.length, burst_data.width, 1)
         create_test_geotiff(burst_data.data_path, dtype='float', value=n, shape=shape)
@@ -38,17 +38,17 @@ class TestMeasurement:
     def test_init(self, burst_datas, gcps):
         measurement = Measurement(burst_datas, gcps, '003.20', 1)
 
-        assert measurement.total_length == 10 * 2
+        assert measurement.total_length == 1000 * 2
         assert measurement.data_mean is None
         assert measurement.data_std is None
 
     def test_get_data(self, burst_datas, gcps):
         measurement = Measurement(burst_datas, gcps, '003.20', 1)
         data = measurement.get_data()
-        assert data.shape == (10 * 2, 20)
+        assert data.shape == (1000 * 2, 2000)
 
-        golden = np.ones((20, 20), dtype=np.complex64)
-        golden[10:, :] *= 2
+        golden = np.ones((2000, 2000), dtype=np.complex64)
+        golden[1000:, :] *= 2
         assert np.allclose(data, golden)
 
     def test_add_metadata(self, burst_datas, gcps):