From 27caa83fbe58626ffd31c2e31ce95022d1e510c0 Mon Sep 17 00:00:00 2001 From: Scott Henderson Date: Thu, 25 Apr 2024 15:30:21 -0700 Subject: [PATCH] convenience scripts for searching ASF --- scripts/findBurstIDs.py | 88 ++++++++++++++++++++++++++++++++ scripts/findSLCs.py | 108 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 196 insertions(+) create mode 100755 scripts/findBurstIDs.py create mode 100755 scripts/findSLCs.py diff --git a/scripts/findBurstIDs.py b/scripts/findBurstIDs.py new file mode 100755 index 0000000..c525abc --- /dev/null +++ b/scripts/findBurstIDs.py @@ -0,0 +1,88 @@ +#! /usr/bin/env python +""" +Search ASF for burstIDs that cover a point + +Usage: findBurstIDs.py LON LAT +Example: findBurstIDs.py -73.604 -49.669 +""" +import argparse +import asf_search as asf +import geopandas as gpd + + +def find_bursts(lon, lat): + """Find ESA Sentinel-1 burstIDs that cover a point""" + results = asf.geo_search( + collections="C2450786986-ASF", + # NOTE: collection doesn't expose beamMode, so need to filter results instead + # beamMode=asf.BEAMMODE.IW, + intersectsWith=f"POINT({lon} {lat})", + ) + gf = gpd.GeoDataFrame.from_features(results.geojson(), crs=4326) + + gf = gf[gf.fileID.str.contains("_IW")].reset_index(drop=True) + # Standard Name (PATH_ID_SWATH) + gf["burstID"] = gf.fileID.str[3:-9] + gf = gf.dropna(axis='columns') + gf = gf.drop(columns='s3Urls') + print(gf.loc[:, ["burstID", "flightDirection"]]) + + return gf + +def slippy_map(gf, lon, lat): + """Plot geopandas polygons using folium and open in browser""" + # Interactive map from CLI + import webbrowser + import folium + from folium.plugins import MiniMap + path = '/tmp/map.html' + m = gf.explore(column='burstID', tiles='Esri.WorldImagery', categorical=True, cmap='Accent') + folium.Marker(location=[lat, lon]).add_to(m) + MiniMap(position="topright", zoom_level_fixed=1).add_to(m) + m.save(path) + print('map saved to', path) + webbrowser.open(f'file://{path}') + + +def static_map(gf, lon, lat): + """Plot geopandas polygons using matplotlib""" + import matplotlib.pyplot as plt + import contextily as cx + + fig, ax = plt.subplots(figsize=(11,8.5)) + gf.plot(ax=ax, column="burstID", facecolor="none", linewidth=4, cmap='Accent', legend=True) + ax.plot(lon, lat, "k*", markersize=10) + + # Zoom out more compared to autoscaling defaults for more context + percent_buffer = .25 + xmin,ymin,xmax,ymax = gf.total_bounds + xbuf = (xmax - xmin) * percent_buffer + ybuf = (ymax - ymin) * percent_buffer + ax.set_xlim(xmin - xbuf, xmax + xbuf) + ax.set_ylim(ymin - ybuf, ymax + ybuf) + # Automatically chosen zoom can be a bit coarse, zoom=9 seems good, or zoom_adjust=1 + cx.add_basemap(ax, source=cx.providers.Esri.WorldImagery, crs='EPSG:4326', zoom_adjust=1) # WorldTerrain + + ax.set_xlabel("Longitude") + ax.set_ylabel("Latitude") + ax.set_title(f"Sentinel-1 Bursts Intersecting ({lon}, {lat})") + plt.save('/tmp/footprints.png') + plt.show() + +if __name__ == "__main__": + # Parse command line arguments + parser = argparse.ArgumentParser( + description="Search ASF for burstIDs that cover a point" + ) + parser.add_argument("lon", type=float, help="Longitude") + parser.add_argument("lat", type=float, help="Latitude") + parser.add_argument("-p", "--show-plot", default=False, action="store_true", help="Plot burstIDs on a map") + #parser.add_argument("-t", "--show-slcs", default=False, action="store_true", help="Lookup all SLCs for burstIDs") + args = parser.parse_args() + #print(args) + gf = find_bursts(args.lon, args.lat) + + if args.show_plot: + #static_map(gf, args.lon, args.lat) + slippy_map(gf, args.lon, args.lat) + diff --git a/scripts/findSLCs.py b/scripts/findSLCs.py new file mode 100755 index 0000000..d6575d1 --- /dev/null +++ b/scripts/findSLCs.py @@ -0,0 +1,108 @@ +#! /usr/bin/env python +""" +Search ASF for SLCs containing a specific burst + +Usage: findSLCs.py BurstID +Example: findSLCs.py 135_289664_IW1 + +""" +import argparse +import asf_search as asf +import geopandas as gpd + + +def get_burst_metadata(burstID): + print('Searching for Burst Metdata...') + results = asf.granule_search([f'S1_{burstID}-BURSTMAP']) + gfB = gpd.GeoDataFrame.from_features(results.geojson(), crs=4326) + gfB["burstID"] = gfB.fileID.str[3:-9] + gfB = gfB.dropna(axis='columns') + print(gfB.iloc[0]) + + return gfB + + +def search_for_slcs(gfB, start, end): + print('Searching for ASF for SLCs...') + burst = gfB.iloc[0] + results = asf.geo_search( + platform=[asf.PLATFORM.SENTINEL1], + processingLevel='SLC', #or BURST from 2023 onwards for select paths + beamMode=asf.BEAMMODE.IW, + relativeOrbit=burst.pathNumber, + intersectsWith=f"POINT({burst.centerLon} {burst.centerLat})", + start=start, + end=end, + ) + + gf = gpd.GeoDataFrame.from_features(results.geojson(), crs=4326) + gf['datetime'] = gpd.pd.to_datetime(gf['startTime']) + + print('BurstID:', burst.burstID) + print('Number of SLCs:', len(gf)) + print('Timespan:', gf.startTime.iloc[0], gf.startTime.iloc[-1]) + print('Polarizations:', list(gf.polarization.unique())) + print('platforms:', list(gf.platform.unique())) + s = gf.datetime.diff(-1).dt.round('1d').dt.days.dropna().astype('i2') + print('Temporal separation (min, mode, max days):', s.min(), s.mode()[0], s.max()) + + return gf + + +def slippy_map(gf, gfB): + """Plot geopandas polygons using folium and open in browser""" + # Interactive map from CLI + import webbrowser + import folium + from folium.plugins import MiniMap + path = '/tmp/mapSLCs.html' + m = gf.drop(columns='datetime').explore(tiles='Esri.WorldImagery', color='cyan', style_kwds=dict(fill=None)) + gfB.explore(m=m, column='burstID', cmap=['magenta'], legend=True) + MiniMap(position="topright").add_to(m) + m.save(path) + print('map saved to', path) + webbrowser.open(f'file://{path}') + + +def timeline(gf, burstID): + """ Plot timeline of SLC acquisitions """ + # Interactive map from CLI + import matplotlib.pyplot as plt + plt.figure(figsize=(11,4)) + plt.scatter(gf.datetime, gf.platform, marker='|', color='k') + plt.title(burstID) + plt.plot() + plt.show() + #plt.savefig(f'/tmp/{burstID}-timeline.pdf') + + +if __name__ == "__main__": + # Parse command line arguments + parser = argparse.ArgumentParser( + description="Search ASF for burstIDs that cover a point" + ) + parser.add_argument("burstID", type=str, help="BurstID (e.g. 106_227373_IW2)") + parser.add_argument("-s", "--start", default=None, type=str, help="Start date (e.g. 2017-01-01)") + parser.add_argument("-e", "--end", default=None, type=str, help="End date (e.g. 2023-01-01)") + parser.add_argument("-g", "--geojson", default=False, action="store_true", help="Save GeoJSON metadata") + parser.add_argument("-p", "--show-plots", default=False, action="store_true", help="Show map of SLC footprints") + + args = parser.parse_args() + # if args.start == 'None': + # args.start = None + print(args) + + gfB = get_burst_metadata(args.burstID) + gf = search_for_slcs(gfB, args.start, args.end) + + # Dump entire list of SLCs + print('------------------------------------') + print('\n'.join(gf.sceneName.to_list())) + + if args.geojson: + # Drop lists before saving to GeoJSON + gf.drop(columns='s3Urls').to_file(f'/tmp/{args.burstID}.geojson', driver='GeoJSON') + + if args.show_plots: + slippy_map(gf, gfB) + timeline(gf, args.burstID)