From e01a6fc017eea6f3151c598b3629cb4b8bb2e238 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 19 Apr 2024 09:04:01 -0500 Subject: [PATCH] correclty compute burst byte offests --- README.md | 7 +++---- src/burst2safe/measurement.py | 24 +++++++++++++++++++++++- src/burst2safe/product.py | 14 ++++++++++++++ src/burst2safe/swath.py | 1 + tests/test_measurement.py | 12 ++++++------ 5 files changed, 47 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index a2e8db8..b26b43b 100644 --- a/README.md +++ b/README.md @@ -67,11 +67,11 @@ 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. @@ -79,7 +79,6 @@ A full accounting of omitted datasets and differing fields can be found below: - 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. diff --git a/src/burst2safe/measurement.py b/src/burst2safe/measurement.py index 2db9d1b..3cc6b5a 100644 --- a/src/burst2safe/measurement.py +++ b/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) diff --git a/src/burst2safe/product.py b/src/burst2safe/product.py index 9da5d66..eea7209 100644 --- a/src/burst2safe/product.py +++ b/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') diff --git a/src/burst2safe/swath.py b/src/burst2safe/swath.py index fe9469a..6f33bfb 100644 --- a/src/burst2safe/swath.py +++ b/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) diff --git a/tests/test_measurement.py b/tests/test_measurement.py index d244893..f32b70f 100644 --- a/tests/test_measurement.py +++ b/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):