initial version of s1_info cli tool (#83)

scottstanie · web-flow · commit 4a6291a98377 · 2022-11-17T18:51:07.000-05:00
* initial version of s1_info cli tool

* print output name

* fix annotation, print path

* add example, account for possible unzipped SAFE dirs

* vary examples to show SAFE dir is okay
diff --git a/setup.cfg b/setup.cfg
@@ -25,3 +25,6 @@ install_requires =
 [options.packages.find]
 where = src
 
+[options.entry_points]
+console_scripts =
+    s1_info = s1reader.s1_info:main
diff --git a/src/s1reader/s1_info.py b/src/s1reader/s1_info.py
@@ -0,0 +1,153 @@
+"""Extract the burst ID information from a Sentinel-1 SLC product."""
+import argparse
+from itertools import chain
+from pathlib import Path
+from typing import List, Optional, Union
+import warnings
+
+import s1reader
+
+
+def get_bursts(
+    filename: Union[Path, str], pol: str = "vv", iw: Optional[int] = None
+) -> List[s1reader.Sentinel1BurstSlc]:
+    if iw is not None:
+        iws = [iw]
+    else:
+        iws = [1, 2, 3]
+    burst_nested_list = [
+        s1reader.load_bursts(filename, None, iw, pol, flag_apply_eap=False)
+        for iw in iws
+    ]
+    return list(chain.from_iterable(burst_nested_list))
+
+
+def _is_safe_dir(path):
+    # Rather than matching the name, we just check for the existence of the
+    # manifest.safe file and annotation files
+    if not (path / "manifest.safe").is_file():
+        return False
+    annotation_dir = path / "annotation"
+    if not annotation_dir.is_dir():
+        return False
+    if len(list(annotation_dir.glob("*.xml"))) == 0:
+        return False
+    return True
+
+
+def _plot_bursts(safe_path, output_dir="burst_maps"):
+    from s1reader.utils import plot_bursts
+
+    orbit_dir = None
+    xs, ys = 5, 10
+    epsg = 4326
+    d = Path(output_dir).resolve()
+    d.mkdir(exist_ok=True, parents=True)
+    print(f"Output directory: {d}")
+    output_filename = d / safe_path.stem
+    plot_bursts.burst_map(safe_path, orbit_dir, xs, ys, epsg, output_filename)
+
+
+EXAMPLE = """
+Example usage:
+
+    # Print all bursts in a Sentinel-1 SLC product
+    s1_info.py S1A_IW_SLC__1SDV_20180601T000000_20180601T000025_021873_025F3D_9E9E.zip
+
+    # Print only the burst IDs
+    s1_info.py S1A_IW_SLC__1SDV_20180601T000000_20180601T000025_021873_025F3D_9E9E.SAFE --burst-id
+
+    # Print burst ids for all files matching the pattern
+    s1_info.py -b S1A_IW_SLC__1SDV_2018*
+
+    # Print only from subswath IW1, and "vv" polarization
+    s1_info.py -b S1A_IW_SLC__1SDV_2018* --iw 1 --pol vv
+
+    # Get info for all products in the 'data/' directory
+    s1_info.py data/
+ 
+    # Plot the burst map, saving files into the 'burst_maps/' directory
+    s1_info.py S1A_IW_SLC__1SDV_20180601T000000_20180601T000025_021873_025F3D_9E9E.SAFE/ --plot
+    s1_info.py S1A_IW_SLC__1SDV_20180601T000000_20180601T000025_021873_025F3D_9E9E.zip -p -o my_burst_maps
+"""
+
+
+def get_cli_args():
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        description="Extract the burst ID information from a Sentinel-1 SLC product.",
+        epilog=EXAMPLE,
+    )
+    parser.add_argument(
+        "paths",
+        help="Path to the Sentinel-1 SLC product(s), or directory containing products.",
+        nargs="+",
+    )
+    parser.add_argument(
+        "--pol",
+        default="vv",
+        choices=["vv", "vh", "hh", "hv"],
+        help="Polarization to use.",
+    )
+    parser.add_argument(
+        "-i",
+        "--iw",
+        type=int,
+        choices=[1, 2, 3],
+        help="Print only the burst IDs for the given IW.",
+    )
+    parser.add_argument(
+        "-b",
+        "--burst-id",
+        action="store_true",
+        help="Print only the burst IDs for all bursts.",
+    )
+    parser.add_argument(
+        "-p",
+        "--plot",
+        action="store_true",
+        help="Plot the burst map for all bursts.",
+    )
+    parser.add_argument(
+        "-o",
+        "--output-dir",
+        default="burst_maps",
+        help=(
+            "Name of the output directory for the burst maps (if plotting),"
+            " with files named for each S1 product (default= %(default)s)."
+        ),
+    )
+    return parser.parse_args()
+
+
+def main():
+    args = get_cli_args()
+    paths = [Path(p) for p in args.paths]
+    all_files = []
+    for path in paths:
+        if path.is_file() or _is_safe_dir(path):
+            all_files.append(path)
+        elif path.is_dir():
+            # Get all matching files within the directory
+            files = path.glob("S1[AB]_IW*")
+            all_files.extend(list(sorted(files)))
+        else:
+            warnings.warn(f"{path} is not a file or directory. Skipping.")
+
+    print(f"Found {len(all_files)} Sentinel-1 SLC products.")
+    for path in all_files:
+        if args.plot:
+            _plot_bursts(path)
+            continue
+        print(f"Bursts in {path}:")
+        print("-" * 80)
+        # Do we want to pretty-print this with rich?
+        for burst in get_bursts(path, args.pol, args.iw):
+            if args.burst_id:
+                print(burst.burst_id)
+            else:
+                print(burst)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/s1reader/utils/plot_bursts.py b/src/s1reader/utils/plot_bursts.py
@@ -1,17 +1,14 @@
 import argparse
 import os
-import sys
-from importlib import import_module
+from pathlib import Path
 
-named_libs = [('fiona', 'fiona'), ('geopandas', 'gpd'), ('pandas', 'pd')]
-
-for (name, short) in named_libs:
-    try:
-        lib = import_module(name)
-    except:
-        print(sys.exc_info())
-    else:
-        globals()[short] = lib
+try:
+    import fiona
+    import geopandas as gpd
+    import pandas as pd
+except ImportError:
+    print("ERROR: fiona, geopandas, and pandas are required for this script.")
+    raise
 
 from osgeo import osr
 from shapely.geometry import Polygon
@@ -25,34 +22,27 @@ def command_line_parser():
     '''
     Command line parser
     '''
-    parser = argparse.ArgumentParser(description="""
-                                     Create a burst map for a single slc""",
-                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
-    parser.add_argument('-s', '--slc', type=str, action='store',
-                        dest='slc',
-                        help="slc to map")
-    parser.add_argument('-d', '--orbit-dir', type=str, dest='orbit_dir',
-                        help="Directory containing orbit files")
+    parser = argparse.ArgumentParser(
+        description="Create a burst map for a single slc",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+    parser.add_argument('-s', '--slc', help="Sentinel-1 product to load.")
+    parser.add_argument('-d', '--orbit-dir', help="Directory containing orbit files")
     parser.add_argument('-x', '--x-spacing', type=float, default=5,
-                        dest='x_spacing',
                         help='Spacing of the geogrid in x direction')
     parser.add_argument('-y', '--y-spacing', type=float, default=10,
-                        dest='y_spacing',
                         help='Spacing of the geogrid in y direction')
-    parser.add_argument('-e', '--epsg', type=int, dest='epsg',
-                        help='EPSG for output coordinates')
-    parser.add_argument('-o', '--output', type=str, default='burst_map.gpkg',
-                        dest='output',
+    parser.add_argument('-e', '--epsg', type=int, help='EPSG for output coordinates')
+    parser.add_argument('-o', '--output', default='burst_map.gpkg', dest='output',
                         help='Base filename for all output burst map products')
     return parser.parse_args()
 
 
-def burst_map(slc, orbit_dir, x_spacing,
-              y_spacing, epsg,
-              output_filename):
-    """
-    Create a CSV of SLC metadata and plot bursts
-    Parameters:
+def burst_map(slc, orbit_dir, x_spacing, y_spacing, epsg, output_filename):
+    """Create a CSV of SLC metadata and plot bursts.
+
+    Parameters
+    ----------
     slc: str
       Path to SLC file
     orbit_dir: str
@@ -65,11 +55,11 @@ def burst_map(slc, orbit_dir, x_spacing,
       EPSG code for the output coodrdinates
     output_filename: str
       Filename used for the output CSV, shp, html, and kml files
-    
-    Returns:
+
+    Returns
+    -------
     output_filename.csv, output_filename.shp, output_filename.html, output_filename.kml
     """
-    
     # Initialize dictionary that will contain all the info for geocoding
     burst_map = {'burst_id':[], 'length': [], 'width': [],
                  'spacing_x': [], 'spacing_y':[], 'min_x': [],
@@ -78,7 +68,7 @@ def burst_map(slc, orbit_dir, x_spacing,
                  'border':[]}
     i_subswath = [1, 2, 3]
     pol = 'vv'
-    orbit_path = get_orbit_file_from_dir(slc, orbit_dir)
+    orbit_path = get_orbit_file_from_dir(slc, orbit_dir) if orbit_dir else None
 
     for subswath in i_subswath:
         ref_bursts = load_bursts(slc, orbit_path, subswath, pol)
@@ -93,6 +83,8 @@ def burst_map(slc, orbit_dir, x_spacing,
             burst_map['first_valid_sample'].append(burst.first_valid_sample)
             burst_map['last_valid_sample'].append(burst.last_valid_sample)
 
+            # TODO: this will ignore the other border for bursts on the antimeridian.
+            # Should probably turn into a MultiPolygon
             poly = burst.border[0]
             # Give some margin to the polygon
             margin = 0.001
@@ -113,6 +105,7 @@ def burst_map(slc, orbit_dir, x_spacing,
                 tgt_x.append(dummy_x)
                 tgt_y.append(dummy_y)
             
+            # TODO: Get the min/max from the burst database
             if epsg == 4326:
                 x_min = x_spacing * (min(tgt_x) / x_spacing)
                 y_min = y_spacing * (min(tgt_y) / y_spacing)
@@ -130,28 +123,29 @@ def burst_map(slc, orbit_dir, x_spacing,
             burst_map['max_x'].append(x_max)
             burst_map['max_y'].append(y_max)
 
+    out_path = Path(output_filename)
+
     # Save generated burst map as csv
     data = pd.DataFrame.from_dict(burst_map)
-    data.to_csv(f'{output_filename}.csv')
+    data.to_csv(out_path.with_suffix('.csv'))
 
     # Create GeoDataFrame to plot bursts on a map
     df = data
     df['border'] = df['border'].apply(wkt.loads)
-    gdf = gpd.GeoDataFrame(df, crs='epsg:4326')
-    gdf = gdf.rename(columns={'border': 'geometry'}).set_geometry('geometry')
-    
+    gdf = gpd.GeoDataFrame(df.rename(columns={'border': 'geometry'}), crs='epsg:4326')
+
+    gdf.to_file(out_path.with_suffix('.gpkg'), driver='GPKG')
     # Save the GeoDataFrame as a shapefile (some people may prefer the format)
-    gdf.to_file(f'{output_filename}.shp')
-    
+    gdf.to_file(out_path.with_suffix('.shp'))
+
     # Save the GeoDataFrame as a kml
-    kml_path = f'{output_filename}.kml'
-    if os.path.isfile(kml_path):
-        os.remove(kml_path)
-    
+    kml_path = out_path.with_suffix('.kml')
+    if kml_path.exists():
+        kml_path.unlink()
+
     fiona.supported_drivers['KML'] = 'rw'
-    gdf.to_file(f'{output_filename}.kml', driver='KML')
-    
-    
+    gdf.to_file(kml_path, driver='KML')
+
     # Plot bursts on an interactive map
     m = gdf.explore(
         column="burst_id", # make choropleth based on "Burst ID" column
@@ -162,8 +156,8 @@ def burst_map(slc, orbit_dir, x_spacing,
         style_kwds=dict(color="black") # use black outline
        )
 
-    m.save(f'{output_filename}.html')
-    
+    m.save(out_path.with_suffix('.html'))
+
 
 if __name__ == '__main__':
     cmd = command_line_parser()
