From 5c185451f7f557cfa3070a6c0dd1e6e98bcb9a99 Mon Sep 17 00:00:00 2001 From: Spicer Bak Date: Mon, 26 Jan 2026 13:00:17 -0500 Subject: [PATCH 1/4] Add getArgusImagery function for FRF tower orthophoto retrieval Implements issue #18 - adds functions to retrieve Argus imagery from the coastalimaging.erdc.dren.mil server: - getArgusImagery(): Synchronous retrieval of Argus orthophotos - Supports image types: timex, var, snap, brightest, darkest - Rounds time to nearest 30-minute interval - Returns dict with image array, time, epochtime, and metadata - threadGetArgusImagery(): Non-blocking background download version - Returns filename immediately while download proceeds in background Closes #18 Co-Authored-By: Claude Opus 4.5 --- murgtools/__init__.py | 3 + murgtools/getdata/__init__.py | 5 +- murgtools/getdata/getDataFRF.py | 149 +++++++++++++++++++++++++++++++ tests/test_getDataFRF.py | 150 ++++++++++++++++++++++++++++++++ 4 files changed, 306 insertions(+), 1 deletion(-) diff --git a/murgtools/__init__.py b/murgtools/__init__.py index 8b1c5a5..d278cc9 100644 --- a/murgtools/__init__.py +++ b/murgtools/__init__.py @@ -2,6 +2,7 @@ from .getdata import getObs, getDataTestBed, gettime, getnc, removeDuplicatesFromDictionary from .getdata import forecastData, getSatelliteImagery, alt_PlotData +from .getdata import getArgusImagery, threadGetArgusImagery __all__ = [ "getObs", @@ -11,5 +12,7 @@ "removeDuplicatesFromDictionary", "forecastData", "getSatelliteImagery", + "getArgusImagery", + "threadGetArgusImagery", "alt_PlotData", ] diff --git a/murgtools/getdata/__init__.py b/murgtools/getdata/__init__.py index 9a1904e..e8f37e9 100644 --- a/murgtools/getdata/__init__.py +++ b/murgtools/getdata/__init__.py @@ -5,7 +5,8 @@ the USACE Field Research Facility (FRF) Coastal Model Test Bed (CMTB). """ -from .getDataFRF import getObs, getDataTestBed, gettime, getnc, removeDuplicatesFromDictionary +from .getDataFRF import (getObs, getDataTestBed, gettime, getnc, removeDuplicatesFromDictionary, + getArgusImagery, threadGetArgusImagery) from .getOutsideData import forecastData, getSatelliteImagery from .getPlotData import alt_PlotData @@ -21,5 +22,7 @@ "removeDuplicatesFromDictionary", "forecastData", "getSatelliteImagery", + "getArgusImagery", + "threadGetArgusImagery", "alt_PlotData", ] diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py index a4663da..224f622 100755 --- a/murgtools/getdata/getDataFRF.py +++ b/murgtools/getdata/getDataFRF.py @@ -3319,3 +3319,152 @@ def getCSHOREOutput(self, prefix): 'runupMean': ncfile['runupMean'][dataIndex], 'runup2perc': ncfile['runup2perc'][dataIndex]} return mod + + +def getArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=True): + """Retrieve Argus orthophoto imagery from the FRF coastal imaging server. + + Argus images are available every 30 minutes from the FRF tower. This function + downloads GeoTIFF orthophotos from the coastalimaging.erdc.dren.mil server. + + Args: + dateOfInterest (datetime): Target datetime for image retrieval. + Will be rounded to nearest 30-minute interval. + filename (str, optional): Path to save the GeoTIFF file. If None, + returns image in memory only. Defaults to None. + imageType (str, optional): Type of Argus image product. Options are: + 'timex' (time exposure average, default), 'var' (variance), + 'snap' (snapshot), 'brightest', 'darkest'. + verbose (bool, optional): Enable logging output. Defaults to True. + + Returns: + dict: Dictionary containing: + - 'image': numpy.ndarray (H, W, 3) uint8 RGB image + - 'time': datetime object of image capture time + - 'epochtime': float, seconds since 1970-01-01 + - 'imageType': str, image type used + - 'filename': str, path if saved (None otherwise) + - 'url': str, source URL + Returns None if image not found. + + Example: + >>> import datetime as DT + >>> result = getArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0)) + >>> if result: + ... print(result['time'], result['image'].shape) + + """ + import requests + import tifffile + import tempfile + + if verbose: + logging.basicConfig(level=logging.INFO) + + # Validate imageType + valid_types = ['timex', 'var', 'snap', 'brightest', 'darkest'] + if imageType.lower() not in valid_types: + raise ValueError(f"Invalid imageType '{imageType}'. Must be one of: {valid_types}") + + # Round to nearest 30 minutes + minutes = dateOfInterest.minute + if minutes < 15: + rounded_minutes = 0 + elif minutes < 45: + rounded_minutes = 30 + else: + rounded_minutes = 0 + dateOfInterest = dateOfInterest + DT.timedelta(hours=1) + + roundedTime = dateOfInterest.replace(minute=rounded_minutes, second=0, microsecond=0) + + # Construct URL + baseURL = "https://coastalimaging.erdc.dren.mil/FrfTower/Processed/Orthophotos/cxgeo/" + fldr = roundedTime.strftime("%Y_%m_%d") + fname = f'{roundedTime.strftime("%Y%m%dT%H%M%SZ")}.FrfTower.cxgeo.{imageType}.tif' + url = urljoin(baseURL, fldr, fname) + + logging.info(f"Retrieving Argus imagery from {url}") + + # Download the image + try: + resp = requests.get(url, stream=True, timeout=60) + resp.raise_for_status() + except requests.exceptions.RequestException as e: + logging.warning(f"Failed to retrieve Argus imagery: {e}") + return None + + # Save to file or temp file + if filename is not None: + with open(filename, 'wb') as f: + for chunk in resp.iter_content(chunk_size=8192): + f.write(chunk) + image = tifffile.imread(filename) + else: + with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp: + tmp_path = tmp.name + for chunk in resp.iter_content(chunk_size=8192): + tmp.write(chunk) + try: + image = tifffile.imread(tmp_path) + finally: + os.unlink(tmp_path) + + # Ensure uint8 RGB format + if image.dtype != np.uint8: + image = np.clip(image / image.max() * 255, 0, 255).astype(np.uint8) + + # Compute epoch time + epochtime = nc.date2num(roundedTime, 'seconds since 1970-01-01') + + logging.info(f"Retrieved Argus {imageType} image: {image.shape}") + + return { + 'image': image, + 'time': roundedTime, + 'epochtime': epochtime, + 'imageType': imageType, + 'filename': filename, + 'url': url, + } + + +def threadGetArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=True): + """Retrieve Argus imagery in a background thread (non-blocking). + + Spawns a daemon thread to download the image and returns immediately + with the output filename. + + Args: + dateOfInterest (datetime): Target datetime for image retrieval. + Will be rounded to nearest 30-minute interval. + filename (str, optional): Path to save the GeoTIFF file. If None, + generates default filename in current directory. + imageType (str, optional): Type of Argus image product. Options are: + 'timex' (default), 'var', 'snap', 'brightest', 'darkest'. + verbose (bool, optional): Enable logging output. Defaults to True. + + Returns: + str: Output filename where image will be saved. + + Example: + >>> import datetime as DT + >>> filename = threadGetArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0)) + >>> print(f"Downloading to: {filename}") + + """ + import threading + + if filename is None: + filename = os.path.join( + os.getcwd(), f'Argus_{imageType}_{dateOfInterest.strftime("%Y%m%dT%H%M%SZ")}.tif' + ) + + t = threading.Thread( + target=getArgusImagery, + args=[dateOfInterest], + kwargs={'filename': filename, 'imageType': imageType, 'verbose': verbose}, + daemon=True + ) + t.start() + return filename diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py index d1256f3..42fb7bd 100644 --- a/tests/test_getDataFRF.py +++ b/tests/test_getDataFRF.py @@ -176,3 +176,153 @@ def test_getdatatestbed_has_data_locations(self, mock_date2num, sample_datetime_ assert hasattr(tb, 'FRFdataloc') assert hasattr(tb, 'chlDataLoc') assert hasattr(tb, 'crunchDataLoc') + + +class TestGetArgusImagery: + """Tests for the getArgusImagery function.""" + + def test_invalid_image_type_raises_error(self): + """Test that invalid imageType raises ValueError.""" + from murgtools.getdata.getDataFRF import getArgusImagery + + with pytest.raises(ValueError, match="Invalid imageType"): + getArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0), imageType="invalid") + + def test_valid_image_types_accepted(self): + """Test that all valid image types are accepted.""" + from murgtools.getdata.getDataFRF import getArgusImagery + import requests + + valid_types = ['timex', 'var', 'snap', 'brightest', 'darkest'] + for img_type in valid_types: + # Patch at the requests module level since it's imported locally + with patch.object(requests, 'get') as mock_get: + mock_get.side_effect = requests.exceptions.RequestException("Network disabled") + # Should not raise ValueError for valid image types + result = getArgusImagery( + DT.datetime(2024, 6, 15, 12, 0, 0), + imageType=img_type, + verbose=False + ) + assert result is None # Network error returns None + + def test_time_rounding_to_30_minutes(self): + """Test that times are rounded to nearest 30 minutes.""" + from murgtools.getdata.getDataFRF import getArgusImagery + import requests + + test_cases = [ + (DT.datetime(2024, 6, 15, 12, 10, 0), "120000"), # rounds down to 12:00 + (DT.datetime(2024, 6, 15, 12, 20, 0), "123000"), # rounds up to 12:30 + (DT.datetime(2024, 6, 15, 12, 40, 0), "123000"), # rounds down to 12:30 + (DT.datetime(2024, 6, 15, 12, 50, 0), "130000"), # rounds up to 13:00 + ] + + for input_time, expected_time_str in test_cases: + with patch.object(requests, 'get') as mock_get: + mock_get.side_effect = requests.exceptions.RequestException("Network disabled") + getArgusImagery(input_time, verbose=False) + + # Verify the URL constructed contains the expected time + call_args = mock_get.call_args + assert call_args is not None, f"requests.get was not called for {input_time}" + url = call_args[0][0] + assert expected_time_str in url, f"Expected {expected_time_str} in URL {url}" + + def test_url_construction(self): + """Test that URL is constructed correctly.""" + from murgtools.getdata.getDataFRF import getArgusImagery + import requests + + with patch.object(requests, 'get') as mock_get: + mock_get.side_effect = requests.exceptions.RequestException("Network disabled") + getArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0), imageType="timex", verbose=False) + + call_args = mock_get.call_args + url = call_args[0][0] + + assert "coastalimaging.erdc.dren.mil" in url + assert "FrfTower/Processed/Orthophotos/cxgeo" in url + assert "2024_06_15" in url + assert "20240615T120000Z" in url + assert "timex.tif" in url + + @pytest.mark.slow + def test_successful_image_retrieval(self): + """Test successful image retrieval (requires network).""" + from murgtools.getdata.getDataFRF import getArgusImagery + + # Use a known date when imagery is likely available + result = getArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0)) + + if result is not None: + assert 'image' in result + assert 'time' in result + assert 'epochtime' in result + assert 'imageType' in result + assert 'url' in result + assert isinstance(result['image'], np.ndarray) + assert result['imageType'] == 'timex' + + def test_returns_none_on_network_error(self): + """Test that function returns None on network error.""" + from murgtools.getdata.getDataFRF import getArgusImagery + import requests + + with patch.object(requests, 'get') as mock_get: + mock_get.side_effect = requests.exceptions.RequestException("Connection failed") + result = getArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0), verbose=False) + + assert result is None + + +class TestThreadGetArgusImagery: + """Tests for the threadGetArgusImagery function.""" + + def test_returns_filename_immediately(self): + """Test that function returns filename immediately without blocking.""" + from murgtools.getdata.getDataFRF import threadGetArgusImagery + import time + + with patch('murgtools.getdata.getDataFRF.getArgusImagery') as mock_get: + # Make the mock slow to ensure we're not waiting + def slow_get(*args, **kwargs): + time.sleep(10) + return None + mock_get.side_effect = slow_get + + start = time.time() + filename = threadGetArgusImagery(DT.datetime(2024, 6, 15, 12, 0, 0)) + elapsed = time.time() - start + + # Should return almost immediately (less than 1 second) + assert elapsed < 1.0 + assert filename is not None + assert filename.endswith('.tif') + + def test_generates_default_filename(self): + """Test that default filename is generated correctly.""" + from murgtools.getdata.getDataFRF import threadGetArgusImagery + + with patch('murgtools.getdata.getDataFRF.getArgusImagery'): + filename = threadGetArgusImagery( + DT.datetime(2024, 6, 15, 12, 0, 0), + imageType="var" + ) + + assert "Argus_var_" in filename + assert "20240615T120000Z" in filename + assert filename.endswith('.tif') + + def test_uses_provided_filename(self): + """Test that provided filename is used.""" + from murgtools.getdata.getDataFRF import threadGetArgusImagery + + with patch('murgtools.getdata.getDataFRF.getArgusImagery'): + custom_filename = "/tmp/my_custom_argus.tif" + filename = threadGetArgusImagery( + DT.datetime(2024, 6, 15, 12, 0, 0), + filename=custom_filename + ) + + assert filename == custom_filename From 5ae33572030eddc0de08810770cabdf13ba04007 Mon Sep 17 00:00:00 2001 From: Spicer Bak Date: Mon, 26 Jan 2026 21:10:17 -0500 Subject: [PATCH 2/4] updated repo with test example case --- examples/test_wave_and_imagery.py | 232 ++++++++++++++++++++++++++++++ murgtools/getdata/__init__.py | 3 +- murgtools/getdata/getDataFRF.py | 37 +++++ 3 files changed, 271 insertions(+), 1 deletion(-) create mode 100644 examples/test_wave_and_imagery.py diff --git a/examples/test_wave_and_imagery.py b/examples/test_wave_and_imagery.py new file mode 100644 index 0000000..358ec68 --- /dev/null +++ b/examples/test_wave_and_imagery.py @@ -0,0 +1,232 @@ +""" +Test script for murgtools wave data and imagery functionality. + +Creates a two-panel figure: +- Top: Wave height time series from all available gauges (last 30 days) +- Bottom: Gauge locations over satellite imagery with Argus imagery overlay (50% opacity) +""" +import datetime as DT +import tempfile +import matplotlib.pyplot as plt +import numpy as np + +from murgtools.getdata import getObs, getSatelliteImagery, getArgusImagery, get_geotiff_extent +from murgtools.utils import geoprocess as gp + + +def main(): + # Define time range: last 30 days with clean datetime boundaries + now = DT.datetime.now() + # End at start of today (midnight), start 30 days before that + d2 = DT.datetime(now.year, now.month, now.day, 0, 0, 0) + d1 = d2 - DT.timedelta(days=30) + + print(f"Fetching wave data from {d1.strftime('%Y-%m-%d %H:%M')} to {d2.strftime('%Y-%m-%d %H:%M')}...") + + # Initialize the observation data retriever + obs = getObs(d1, d2) + + # All available wave gauges + gauge_list = [ + 'waverider-26m', + 'waverider-17m', + 'awac-11m', + '8m-array', + 'awac-6m', + 'awac-4.5m', + 'adop-3.5m', + 'xp200m', + 'xp150m', + 'xp125m', + ] + + # Retrieve wave data from all gauges + wave_data = {} + gauge_locations = [] + + for gauge in gauge_list: + print(f" Fetching {gauge}...") + try: + data = obs.getWaveData(gaugenumber=gauge, roundto=30) + if data is not None and 'Hs' in data: + wave_data[gauge] = data + # Get lat/lon from xFRF/yFRF using geoprocess + xFRF = data.get('xFRF') + yFRF = data.get('yFRF') + if xFRF is not None and yFRF is not None: + coords = gp.FRFcoord(xFRF, yFRF, coordType='FRF') + lat = coords['Lat'] + lon = coords['Lon'] + else: + lat = data.get('lat') + lon = data.get('lon') + + gauge_locations.append({ + 'name': gauge, + 'lat': lat, + 'lon': lon, + 'xFRF': xFRF, + 'yFRF': yFRF, + 'depth': data.get('depth'), + }) + print(f" Got {len(data['Hs'])} records") + else: + print(f" No data available") + except Exception as e: + print(f" Error: {e}") + + # Define corners for satellite imagery (wider FRF area to include all gauges) + # Satellite extent: ~2km offshore, 2km alongshore + corners = [ + (36.195, -75.760), # NW - further offshore and north + (36.195, -75.720), # NE + (36.175, -75.760), # SW + (36.175, -75.720), # SE + ] + + # Get satellite imagery (finds closest available to current time) + print("\nFetching satellite imagery...") + sat_result = getSatelliteImagery( + corners, + date=now, + max_cloud_cover=30, + collection='sentinel-2-l2a' + ) + if sat_result: + print(f" Got satellite image from {sat_result['time'].strftime('%Y-%m-%d')}") + print(f" Shape: {sat_result['image'].shape}, Cloud cover: {sat_result['cloud_cover']:.1f}%") + else: + print(" No satellite imagery available") + + # Get Argus imagery (finds closest 30-min interval to current time) + print("\nFetching Argus imagery...") + argus_result = None + argus_extent = None + + # Save to temp file so we can extract extent + with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp: + tmp_path = tmp.name + + try: + argus_result = getArgusImagery(now, filename=tmp_path, imageType='timex') + if argus_result: + print(f" Got Argus image from {argus_result['time'].strftime('%Y-%m-%d %H:%M')}") + print(f" Shape: {argus_result['image'].shape}") + # Extract extent from the GeoTIFF (in State Plane coordinates) + sp_extent = get_geotiff_extent(tmp_path) + print(f" State Plane extent: {sp_extent}") + + # Convert State Plane corners to lat/lon + # sp_extent is [left, right, bottom, top] in State Plane Easting/Northing + sp_left, sp_right, sp_bottom, sp_top = sp_extent + + # Convert corners from State Plane to lat/lon + ll_corner = gp.FRFcoord(sp_left, sp_bottom, coordType='ncsp') # SW corner + ur_corner = gp.FRFcoord(sp_right, sp_top, coordType='ncsp') # NE corner + + # Build lat/lon extent [left, right, bottom, top] as [lon_min, lon_max, lat_min, lat_max] + argus_extent = [ll_corner['Lon'], ur_corner['Lon'], ll_corner['Lat'], ur_corner['Lat']] + print(f" Lat/Lon extent: {argus_extent}") + except Exception as e: + print(f" Error fetching Argus imagery: {e}") + + # Clean up temp file + import os + if os.path.exists(tmp_path): + os.unlink(tmp_path) + + # Create the figure with height ratio - top plot shorter, bottom plot taller + fig, axes = plt.subplots(2, 1, figsize=(14, 12), gridspec_kw={'height_ratios': [1, 2.5]}) + + # ========================================= + # TOP PLOT: Wave height time series + # ========================================= + ax1 = axes[0] + + colors = plt.cm.tab10(np.linspace(0, 1, len(wave_data))) + + for (gauge_name, data), color in zip(wave_data.items(), colors): + if 'time' in data and 'Hs' in data: + ax1.plot(data['time'], data['Hs'], label=gauge_name, color=color, linewidth=1.5) + + ax1.set_xlabel('Time') + ax1.set_ylabel('Significant Wave Height (m)') + ax1.set_title(f'Wave Height at FRF Gauges\n{d1.strftime("%Y-%m-%d")} to {d2.strftime("%Y-%m-%d")}') + ax1.legend(loc='upper right', fontsize=8, ncol=2) + ax1.grid(True, alpha=0.3) + ax1.tick_params(axis='x', rotation=45) + + # ========================================= + # BOTTOM PLOT: Gauge locations over imagery + # ========================================= + ax2 = axes[1] + + # Plot satellite imagery as base layer (larger area) + if sat_result is not None: + ax2.imshow(sat_result['image'], extent=sat_result['extent'], aspect='auto', zorder=1) + + # Overlay Argus imagery at 50% opacity (smaller nearshore inset) + if argus_result is not None and argus_extent is not None: + ax2.imshow(argus_result['image'], extent=argus_extent, aspect='auto', alpha=0.5, zorder=2) + + # Plot gauge locations + for i, gauge in enumerate(gauge_locations): + if gauge['lat'] is not None and gauge['lon'] is not None: + color_idx = gauge_list.index(gauge['name']) % len(colors) + ax2.scatter(gauge['lon'], gauge['lat'], s=100, c=[colors[color_idx]], + edgecolors='white', linewidth=2, zorder=10) + ax2.annotate(gauge['name'], (gauge['lon'], gauge['lat']), + xytext=(5, 5), textcoords='offset points', + fontsize=7, color='white', fontweight='bold', + bbox=dict(boxstyle='round,pad=0.2', facecolor='black', alpha=0.7), + zorder=11) + + ax2.set_xlabel('Longitude') + ax2.set_ylabel('Latitude') + + # Build title based on what imagery is available + title_parts = ['FRF Wave Gauge Locations'] + if sat_result: + title_parts.append(f'Satellite: {sat_result["time"].strftime("%Y-%m-%d")}') + if argus_result: + title_parts.append(f'Argus overlay (50% opacity): {argus_result["time"].strftime("%Y-%m-%d %H:%M")}') + ax2.set_title('\n'.join(title_parts)) + + # Set extent to satellite extent or compute from gauge locations + if sat_result is not None: + ax2.set_xlim(sat_result['extent'][0], sat_result['extent'][1]) + ax2.set_ylim(sat_result['extent'][2], sat_result['extent'][3]) + else: + # Fallback: set extent based on gauge locations + lons = [g['lon'] for g in gauge_locations if g['lon'] is not None] + lats = [g['lat'] for g in gauge_locations if g['lat'] is not None] + if lons and lats: + pad = 0.01 + ax2.set_xlim(min(lons) - pad, max(lons) + pad) + ax2.set_ylim(min(lats) - pad, max(lats) + pad) + + plt.tight_layout() + + # Save the figure + output_file = 'frf_wave_and_imagery.png' + plt.savefig(output_file, dpi=150, bbox_inches='tight') + print(f"\nFigure saved to: {output_file}") + + plt.show() + + # Print summary of gauge locations + print("\n" + "="*80) + print("Gauge Location Summary:") + print("="*80) + for gauge in gauge_locations: + depth_str = f"{gauge['depth']:.1f}m" if gauge['depth'] is not None else "N/A" + xfrf_str = f"{gauge['xFRF']:7.1f}m" if gauge['xFRF'] is not None else "N/A" + yfrf_str = f"{gauge['yFRF']:7.1f}m" if gauge['yFRF'] is not None else "N/A" + lat_str = f"{gauge['lat']:.6f}" if gauge['lat'] is not None else "N/A" + lon_str = f"{gauge['lon']:.6f}" if gauge['lon'] is not None else "N/A" + print(f" {gauge['name']:15s} | lat: {lat_str}, lon: {lon_str} | " + f"xFRF: {xfrf_str}, yFRF: {yfrf_str} | depth: {depth_str}") + + +if __name__ == '__main__': + main() diff --git a/murgtools/getdata/__init__.py b/murgtools/getdata/__init__.py index e8f37e9..8f4445d 100644 --- a/murgtools/getdata/__init__.py +++ b/murgtools/getdata/__init__.py @@ -6,7 +6,7 @@ """ from .getDataFRF import (getObs, getDataTestBed, gettime, getnc, removeDuplicatesFromDictionary, - getArgusImagery, threadGetArgusImagery) + get_geotiff_extent, getArgusImagery, threadGetArgusImagery) from .getOutsideData import forecastData, getSatelliteImagery from .getPlotData import alt_PlotData @@ -22,6 +22,7 @@ "removeDuplicatesFromDictionary", "forecastData", "getSatelliteImagery", + "get_geotiff_extent", "getArgusImagery", "threadGetArgusImagery", "alt_PlotData", diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py index 224f622..236070a 100755 --- a/murgtools/getdata/getDataFRF.py +++ b/murgtools/getdata/getDataFRF.py @@ -3321,6 +3321,43 @@ def getCSHOREOutput(self, prefix): return mod +def get_geotiff_extent(filepath): + """Extract matplotlib extent from GeoTIFF using tifffile. + + Parses GeoTIFF tags (ModelTiepointTag and ModelPixelScaleTag) to compute + geographic bounds suitable for matplotlib imshow extent parameter. + + Args: + filepath (str): Path to GeoTIFF file. + + Returns: + list: [left, right, bottom, top] extent in geographic coordinates + (typically lon/lat or projected coordinates depending on the GeoTIFF). + + Raises: + KeyError: If required GeoTIFF tags are not present in the file. + + Example: + >>> extent = get_geotiff_extent('/path/to/image.tif') + >>> plt.imshow(image, extent=extent) + + """ + import tifffile + + with tifffile.TiffFile(filepath) as tif: + tags = tif.pages[0].tags + # GeoTIFF tags: 33922=ModelTiepointTag, 33550=ModelPixelScaleTag + tiepoint = tags[33922].value # (i, j, k, x, y, z) + scale = tags[33550].value # (scaleX, scaleY, scaleZ) + height, width = tif.pages[0].shape[:2] + # Compute extent: [left, right, bottom, top] + left = tiepoint[3] + top = tiepoint[4] + right = left + width * scale[0] + bottom = top - height * scale[1] + return [left, right, bottom, top] + + def getArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=True): """Retrieve Argus orthophoto imagery from the FRF coastal imaging server. From ade6d05136cbcc0fdff22b6c45dbfe810a8410e1 Mon Sep 17 00:00:00 2001 From: Spicer Bak Date: Thu, 29 Jan 2026 17:04:57 -0500 Subject: [PATCH 3/4] Add findArgusImagery wrapper with search functionality - Add findArgusImagery() wrapper function that searches for nearest available imagery when exact time is not available - Supports two search modes: - method=0: Nearest in TIME (bidirectional search) - method=1: Nearest in HISTORY (backward only, for operational use) - Fix image type names: 'brightest'->'bright', 'darkest'->'dark' - Add search_window_hours parameter (default 24h) - Returns time_requested and time_offset_minutes in result dict - Update test_wave_and_imagery.py to use new wrapper Closes #18 Co-Authored-By: Claude Opus 4.5 --- examples/test_wave_and_imagery.py | 19 ++++-- murgtools/getdata/__init__.py | 3 +- murgtools/getdata/getDataFRF.py | 109 +++++++++++++++++++++++++++++- 3 files changed, 123 insertions(+), 8 deletions(-) diff --git a/examples/test_wave_and_imagery.py b/examples/test_wave_and_imagery.py index 358ec68..1e70492 100644 --- a/examples/test_wave_and_imagery.py +++ b/examples/test_wave_and_imagery.py @@ -10,7 +10,7 @@ import matplotlib.pyplot as plt import numpy as np -from murgtools.getdata import getObs, getSatelliteImagery, getArgusImagery, get_geotiff_extent +from murgtools.getdata import getObs, getSatelliteImagery, findArgusImagery, get_geotiff_extent from murgtools.utils import geoprocess as gp @@ -98,7 +98,7 @@ def main(): else: print(" No satellite imagery available") - # Get Argus imagery (finds closest 30-min interval to current time) + # Get Argus imagery (searches for nearest available within 24-hour window) print("\nFetching Argus imagery...") argus_result = None argus_extent = None @@ -108,10 +108,17 @@ def main(): tmp_path = tmp.name try: - argus_result = getArgusImagery(now, filename=tmp_path, imageType='timex') + # Use findArgusImagery with method=1 (nearest in history) to get most recent image + # Using 'bright' (brightest pixels composite) for best visualization + argus_result = findArgusImagery(now, filename=tmp_path, imageType='bright', + search_window_hours=48, method=1) if argus_result: print(f" Got Argus image from {argus_result['time'].strftime('%Y-%m-%d %H:%M')}") print(f" Shape: {argus_result['image'].shape}") + # Show time offset if search was needed + offset_mins = argus_result.get('time_offset_minutes', 0) + if offset_mins != 0: + print(f" Time offset from requested: {offset_mins} minutes ({offset_mins/60:.1f} hours)") # Extract extent from the GeoTIFF (in State Plane coordinates) sp_extent = get_geotiff_extent(tmp_path) print(f" State Plane extent: {sp_extent}") @@ -189,7 +196,11 @@ def main(): if sat_result: title_parts.append(f'Satellite: {sat_result["time"].strftime("%Y-%m-%d")}') if argus_result: - title_parts.append(f'Argus overlay (50% opacity): {argus_result["time"].strftime("%Y-%m-%d %H:%M")}') + argus_title = f'Argus overlay (50% opacity): {argus_result["time"].strftime("%Y-%m-%d %H:%M")}' + offset_mins = argus_result.get('time_offset_minutes', 0) + if offset_mins != 0: + argus_title += f' (offset: {offset_mins/60:.1f}h)' + title_parts.append(argus_title) ax2.set_title('\n'.join(title_parts)) # Set extent to satellite extent or compute from gauge locations diff --git a/murgtools/getdata/__init__.py b/murgtools/getdata/__init__.py index 8f4445d..9da3d44 100644 --- a/murgtools/getdata/__init__.py +++ b/murgtools/getdata/__init__.py @@ -6,7 +6,7 @@ """ from .getDataFRF import (getObs, getDataTestBed, gettime, getnc, removeDuplicatesFromDictionary, - get_geotiff_extent, getArgusImagery, threadGetArgusImagery) + get_geotiff_extent, getArgusImagery, threadGetArgusImagery, findArgusImagery) from .getOutsideData import forecastData, getSatelliteImagery from .getPlotData import alt_PlotData @@ -25,5 +25,6 @@ "get_geotiff_extent", "getArgusImagery", "threadGetArgusImagery", + "findArgusImagery", "alt_PlotData", ] diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py index 236070a..1953178 100755 --- a/murgtools/getdata/getDataFRF.py +++ b/murgtools/getdata/getDataFRF.py @@ -3371,7 +3371,7 @@ def getArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=Tr returns image in memory only. Defaults to None. imageType (str, optional): Type of Argus image product. Options are: 'timex' (time exposure average, default), 'var' (variance), - 'snap' (snapshot), 'brightest', 'darkest'. + 'snap' (snapshot), 'bright' (brightest pixels), 'dark' (darkest pixels). verbose (bool, optional): Enable logging output. Defaults to True. Returns: @@ -3399,7 +3399,7 @@ def getArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=Tr logging.basicConfig(level=logging.INFO) # Validate imageType - valid_types = ['timex', 'var', 'snap', 'brightest', 'darkest'] + valid_types = ['timex', 'var', 'snap', 'bright', 'dark'] if imageType.lower() not in valid_types: raise ValueError(f"Invalid imageType '{imageType}'. Must be one of: {valid_types}") @@ -3478,7 +3478,7 @@ def threadGetArgusImagery(dateOfInterest, filename=None, imageType="timex", verb filename (str, optional): Path to save the GeoTIFF file. If None, generates default filename in current directory. imageType (str, optional): Type of Argus image product. Options are: - 'timex' (default), 'var', 'snap', 'brightest', 'darkest'. + 'timex' (default), 'var', 'snap', 'bright', 'dark'. verbose (bool, optional): Enable logging output. Defaults to True. Returns: @@ -3505,3 +3505,106 @@ def threadGetArgusImagery(dateOfInterest, filename=None, imageType="timex", verb ) t.start() return filename + + +def findArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=True, + search_window_hours=24, method=0): + """Find available Argus imagery with configurable search strategy. + + Wrapper around getArgusImagery that searches for available imagery + when the exact requested time is not available. + + Args: + dateOfInterest (datetime): Target datetime for image retrieval. + filename (str, optional): Path to save the GeoTIFF file. If None, + returns image in memory only. Defaults to None. + imageType (str, optional): Type of Argus image product. Options are: + 'timex' (time exposure average, default), 'var' (variance), + 'snap' (snapshot), 'bright' (brightest pixels), 'dark' (darkest pixels). + verbose (bool, optional): Enable logging output. Defaults to True. + search_window_hours (int, optional): Hours to search for available + imagery. Default 24 hours. + method (int, optional): Search strategy. Options are: + 0 = Nearest in TIME (bidirectional, returns closest available) + 1 = Nearest in HISTORY (backward only, returns most recent before target) + Defaults to 0. + + Returns: + dict: Same as getArgusImagery, plus: + - 'time_requested': original requested datetime (rounded to 30-min) + - 'time_offset_minutes': offset from requested time (negative=earlier) + Returns None if no image found within search window. + + Example: + >>> import datetime as DT + >>> # Find nearest available image to now + >>> result = findArgusImagery(DT.datetime.now(), method=0) + >>> if result: + ... print(f"Found image {result['time_offset_minutes']} min from requested") + >>> # Find most recent image before target (for operational use) + >>> result = findArgusImagery(DT.datetime.now(), method=1) + + """ + if verbose: + logging.basicConfig(level=logging.INFO) + + # Round to nearest 30 minutes (same logic as getArgusImagery) + minutes = dateOfInterest.minute + if minutes < 15: + rounded_minutes = 0 + rounded_time = dateOfInterest + elif minutes < 45: + rounded_minutes = 30 + rounded_time = dateOfInterest + else: + rounded_minutes = 0 + rounded_time = dateOfInterest + DT.timedelta(hours=1) + time_requested = rounded_time.replace(minute=rounded_minutes, second=0, microsecond=0) + + # Calculate max number of 30-min slots to search + max_slots = int(search_window_hours * 2) + + # Generate candidate timestamps based on method + candidates = [time_requested] # Always try exact time first + + if method == 0: + # Nearest in TIME: alternating offsets [-30, +30, -60, +60, ...] + for i in range(1, max_slots + 1): + offset_minutes = i * 30 + candidates.append(time_requested - DT.timedelta(minutes=offset_minutes)) + candidates.append(time_requested + DT.timedelta(minutes=offset_minutes)) + elif method == 1: + # Nearest in HISTORY: backward only [-30, -60, -90, ...] + for i in range(1, max_slots + 1): + offset_minutes = i * 30 + candidates.append(time_requested - DT.timedelta(minutes=offset_minutes)) + else: + raise ValueError(f"Invalid method {method}. Must be 0 (nearest in time) or 1 (nearest in history).") + + # Try each candidate until we find one + for candidate_time in candidates: + if verbose: + logging.info(f"Trying Argus imagery for {candidate_time.strftime('%Y-%m-%d %H:%M')}...") + + result = getArgusImagery(candidate_time, filename=filename, imageType=imageType, verbose=False) + + if result is not None: + # Calculate offset from requested time + time_offset = candidate_time - time_requested + time_offset_minutes = int(time_offset.total_seconds() / 60) + + # Add search metadata to result + result['time_requested'] = time_requested + result['time_offset_minutes'] = time_offset_minutes + + if verbose: + logging.info(f"Found Argus {imageType} image at {candidate_time.strftime('%Y-%m-%d %H:%M')} " + f"(offset: {time_offset_minutes} minutes)") + + return result + + # No image found within search window + if verbose: + logging.warning(f"No Argus imagery found within {search_window_hours} hours of " + f"{time_requested.strftime('%Y-%m-%d %H:%M')}") + return None From 9e1452ce879b0bf93be409fa550f161245303795 Mon Sep 17 00:00:00 2001 From: Spicer Bak Date: Thu, 29 Jan 2026 17:06:54 -0500 Subject: [PATCH 4/4] Fix test to use correct Argus image type names Update valid_types in test to match actual server naming: - 'brightest' -> 'bright' - 'darkest' -> 'dark' Co-Authored-By: Claude Opus 4.5 --- tests/test_getDataFRF.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py index 42fb7bd..e121f99 100644 --- a/tests/test_getDataFRF.py +++ b/tests/test_getDataFRF.py @@ -193,7 +193,7 @@ def test_valid_image_types_accepted(self): from murgtools.getdata.getDataFRF import getArgusImagery import requests - valid_types = ['timex', 'var', 'snap', 'brightest', 'darkest'] + valid_types = ['timex', 'var', 'snap', 'bright', 'dark'] for img_type in valid_types: # Patch at the requests module level since it's imported locally with patch.object(requests, 'get') as mock_get: