From 4323a90cc2e0cf6263f30dd08af8a5d53bc7a9f2 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Mon, 27 Jul 2026 15:30:50 -0700 Subject: [PATCH 01/21] gracefully handle missing availability service --- PIQQA.py | 230 ++++++++------ reportUtils.py | 801 ++++++++++++++++++++++++++++++------------------- 2 files changed, 635 insertions(+), 396 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index cafbb9b..6e06f10 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -39,6 +39,7 @@ import numpy as np import shutil import math +import urllib.request def checkAvailability( @@ -125,7 +126,29 @@ def doAvailability( nBottom, tolerance, imageDir, + tsMetricsExist, ): + # test out whether the availability service is up and running, if not, then skip the availability plot + availability_exists = True + availablity_test_url = "https://service.earthscope.org/fdsnws/availability/1" + # if the service is returning an HTTP 410, then bypass the availablity plot + try: + urllib.request.urlopen(availablity_test_url) + except urllib.error.HTTPError as e: + if e.code == 410: + print("WARNING: Availability service is down, skipping availability plot") + else: + print( + f"WARNING: Availability service returned HTTP {e.code}, skipping availability plot" + ) + availability_exists = False + return {}, splitPlots, {}, [], "", availability_exists + except Exception as e: + print( + f"WARNING: Availability service is down ({e}), skipping availability plot" + ) + availability_exists = False + return {}, splitPlots, {}, [], "", availability_exists print("INFO: Producing availability plot") avFilesDictionary = {} # collects plot filenames @@ -162,30 +185,32 @@ def doAvailability( tmpMetadataDF = allMetadataDF[allMetadataDF["channel"] == channelGroup] ## Get percent availability numbers for all targets for this channelGroup so that you can narrow it down to the top/bottom available - # First try to use ts_percent_availability_total + # First try to use ts_percent_availability_total, unless we already know ts_ metrics are down pctAvDF = pd.DataFrame() - print( - f" INFO: Retrieving percent availability information for {channelGroup}" - ) - try: - pctAvDF = reportUtils.addMetricToDF( - "ts_percent_availability_total", - pctAvDF, - network, - stations, - locations, - channelGroup, - startDate, - endDate, + if tsMetricsExist: + print( + f" INFO: Retrieving percent availability information for {channelGroup}" ) - except: - print(f"ERROR: Trouble accessing ts_percent_availability_total") + try: + pctAvDF = reportUtils.addMetricToDF( + "ts_percent_availability_total", + pctAvDF, + network, + stations, + locations, + channelGroup, + startDate, + endDate, + ) + except: + print(f"ERROR: Trouble accessing ts_percent_availability_total") - if not pctAvDF.empty: - pctAvDF.rename( - columns={"ts_percent_availability_total": "availability"}, inplace=True - ) + if not pctAvDF.empty: + pctAvDF.rename( + columns={"ts_percent_availability_total": "availability"}, + inplace=True, + ) # If not all targets have ts_percent_availability_total, then try getting it from percent_availability @@ -405,6 +430,7 @@ def doAvailability( bottomStationsAvDF, services = reportUtils.getAvailability( bottomStations["snclq"], startDate, endDate, tolerance, "" ) + for service in services: if service not in services_used: services_used.append(service) @@ -422,7 +448,6 @@ def doAvailability( stn = 0 for station in bottomStationList: - # data extents and gaps extents if ( not bottomStations[ @@ -715,6 +740,7 @@ def doAvailability( avFilesDictionary, services_used, availabilityType, + availability_exists, ) @@ -736,10 +762,17 @@ def doBoxPlots( nBottom, includeOutliers, imageDir, + tsMetricsExist, ): # Retrieve metrics from mustang service print("INFO: Generating Boxplots") print(" INFO: Retrieving metrics from the MUSTANG webservices...") + if not tsMetricsExist: + if "ts_num_gaps_total" in metricList: + metricList.remove("ts_num_gaps_total") + if "ts_num_gaps" in metricList: + metricList.remove("ts_num_gaps") + metricList.append("num_gaps") metricsWithPlots = list() metricDF = pd.DataFrame() @@ -1140,6 +1173,7 @@ def doPDFs( network, stations, locations, + availabilityExists, ): pdfDictionary = {} @@ -1194,18 +1228,17 @@ def doPDFs( thisChannel = lowestTarget.split(".")[3] # Check the availablity service to get the start and end to the data - thisStationAvDF, service = reportUtils.getAvailability( - [lowestTarget], startDate, endDate, tolerance, "" - ) - thisAvail = thisStationAvDF[ - (thisStationAvDF["station"] == thisStation) - & (thisStationAvDF["network"] == thisNetwork) - & (thisStationAvDF["location"] == thisLocation) - & (thisStationAvDF["channel"] == thisChannel) - ] - - thisStart = thisAvail["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() - thisEnd = thisAvail["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() + thisStart, thisEnd = reportUtils.getDataStartEnd( + lowestTarget, + thisNetwork, + thisStation, + thisLocation, + thisChannel, + startDate, + endDate, + tolerance, + availabilityExists, + ) print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" @@ -1268,18 +1301,17 @@ def doPDFs( thisChannel = largestTarget.split(".")[3] # Get start and end to the data for this target - thisStationAvDF, service = reportUtils.getAvailability( - [largestTarget], startDate, endDate, tolerance, "" - ) - thisAvail = thisStationAvDF[ - (thisStationAvDF["station"] == thisStation) - & (thisStationAvDF["network"] == thisNetwork) - & (thisStationAvDF["location"] == thisLocation) - & (thisStationAvDF["channel"] == thisChannel) - ] - - thisStart = thisAvail["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() - thisEnd = thisAvail["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() + thisStart, thisEnd = reportUtils.getDataStartEnd( + largestTarget, + thisNetwork, + thisStation, + thisLocation, + thisChannel, + startDate, + endDate, + tolerance, + availabilityExists, + ) print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" @@ -1378,6 +1410,7 @@ def doSpectrograms( powerRanges, spectColorPalette, imageDir, + availabilityExists, ): spectDictionary = {} # used to track the plots to be used in the final report @@ -1444,18 +1477,17 @@ def doSpectrograms( thisChannel = lowestTarget.split(".")[3] # Get the start and end to the data for this target - thisStationAvDF, service = reportUtils.getAvailability( - [lowestTarget], startDate, endDate, tolerance, "" - ) - thisAvail = thisStationAvDF[ - (thisStationAvDF["station"] == thisStation) - & (thisStationAvDF["network"] == thisNetwork) - & (thisStationAvDF["location"] == thisLocation) - & (thisStationAvDF["channel"] == thisChannel) - ] - - thisStart = thisAvail["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() - thisEnd = thisAvail["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() + thisStart, thisEnd = reportUtils.getDataStartEnd( + lowestTarget, + thisNetwork, + thisStation, + thisLocation, + thisChannel, + startDate, + endDate, + tolerance, + availabilityExists, + ) print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" @@ -1517,18 +1549,17 @@ def doSpectrograms( thisChannel = largestTarget.split(".")[3] # Get the start/end times of the data for this target - thisStationAvDF, service = reportUtils.getAvailability( - [largestTarget], startDate, endDate, tolerance, "" - ) - thisAvail = thisStationAvDF[ - (thisStationAvDF["station"] == thisStation) - & (thisStationAvDF["network"] == thisNetwork) - & (thisStationAvDF["location"] == thisLocation) - & (thisStationAvDF["channel"] == thisChannel) - ] - - thisStart = thisAvail["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() - thisEnd = thisAvail["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() + thisStart, thisEnd = reportUtils.getDataStartEnd( + largestTarget, + thisNetwork, + thisStation, + thisLocation, + thisChannel, + startDate, + endDate, + tolerance, + availabilityExists, + ) print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" @@ -1704,6 +1735,7 @@ def doReport( spectColorPalette, powerRanges, userDefinedPowerRange, + availabilityExists, ): print(f"INFO: Writing to file {outfile}") @@ -1749,11 +1781,17 @@ def doReport( ) ###### AVAILABILITY ###### - for service in services: - stationServiceLink = f'{service}-station service,' - availabilityServiceLink = f'{service}-availability service,' - stationServiceLink = f"{stationServiceLink[:-1]}" - availabilityServiceLink = f"{availabilityServiceLink[:-1]}" + # for service in services: + # stationServiceLink = f'{service}-station service,' + # availabilityServiceLink = f'{service}-availability service,' + # stationServiceLink = f"{stationServiceLink[:-1]}" + if availabilityExists: + availabilityServiceLink = ( + f"http://service.earthscope.org/fdsnws/availability/1/" + ) + else: + availabilityServiceLink = f'mustang service' + stationServiceLink = f'fdsnws-station service' availabilityIntroText = ( f"The plot(s) below gives an overview of the available data for the requested timespan. " @@ -2753,22 +2791,32 @@ def main(): nTop = int(nBoxPlotSta / 2) nBottom = int(nBoxPlotSta - nTop) + # Check independently whether the ts_ (time-series) MUSTANG metrics exist right now - + # this can go down (or come back) on its own schedule, separate from the fdsnws-availability service + tsMetricsExist = reportUtils.checkTsMetricsExist(startDate, endDate) + # Create Availability Plots - [topStationsDict, splitPlots, avFilesDictionary, services, availabilityType] = ( - doAvailability( - splitPlots, - startDate, - endDate, - network, - stations, - locations, - channels, - nBoxPlotSta, - nTop, - nBottom, - tolerance, - imageDir, - ) + [ + topStationsDict, + splitPlots, + avFilesDictionary, + services, + availabilityType, + availabilityExists, + ] = doAvailability( + splitPlots, + startDate, + endDate, + network, + stations, + locations, + channels, + nBoxPlotSta, + nTop, + nBottom, + tolerance, + imageDir, + tsMetricsExist, ) # Create BoxPlots @@ -2788,6 +2836,7 @@ def main(): nBottom, includeOutliers, imageDir, + tsMetricsExist, ) ) @@ -2806,6 +2855,7 @@ def main(): network, stations, locations, + availabilityExists, ) # Create Spectrogram Plots @@ -2821,6 +2871,7 @@ def main(): powerRanges, spectColorPalette, imageDir, + availabilityExists, ) # Generate Station Map @@ -2868,6 +2919,7 @@ def main(): spectColorPalette, powerRanges, userDefinedPowerRange, + availabilityExists, ) # Zip the report diff --git a/reportUtils.py b/reportUtils.py index 9406d3a..0b8449a 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -1,8 +1,8 @@ -''' +""" :copyright: Copyright (C) 2021 Laura Keyson, IRIS Data Management Center - -:license: + +:license: Licensed to the Apache Software Foundation (ASF) under one or more contributor license agreements. See the NOTICE file distributed with this work for additional information @@ -18,10 +18,9 @@ "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations -under the License. - -''' +under the License. +""" import requests from io import StringIO @@ -31,451 +30,639 @@ import urllib +def get_extents_from_mustang(snclqs, startDate, endDate): + availabilityDF = pd.DataFrame() + for snclq in snclqs: + staloc = ( + snclq.split(".")[1] + + "." + + (snclq.split(".")[2] if snclq.split(".")[2] != "" else "--") + ) + URL = ( + f"https://service.earthscope.org/mustang/measurements/1/" + f"query?metric=percent_availability&target={snclq}&format=text&timewindow={startDate},{endDate}&nodata=404" + ) + + try: + tempDF = pd.read_csv(URL, header=1, parse_dates=["start", "end"]) + + # find the first and last non-0 availability values and use those as the start and end times + nonZeroDF = tempDF[tempDF["value"] > 0] + if not nonZeroDF.empty: + startTime = nonZeroDF["start"].min() + endTime = nonZeroDF["end"].max() + tempDF = pd.DataFrame( + { + "network": [snclq.split(".")[0]], + "station": [snclq.split(".")[1]], + "location": [snclq.split(".")[2]], + "channel": [snclq.split(".")[3]], + "quality": [snclq.split(".")[4]], + "staloc": [staloc], + "earliest": [startTime], + "latest": [endTime], + } + ) + else: + tempDF = pd.DataFrame() + + availabilityDF = pd.concat([availabilityDF, tempDF], ignore_index=True) + except Exception as e: + print( + f" INFO: Unable to retrieve availability for {snclq} from Mustang service" + ) + print(f" {URL}") + pass + + return availabilityDF + + def getAvailability(snclqs, startDate, endDate, tolerance, avtype): availabilityDF = pd.DataFrame() services = [] - - - if avtype == '': + + if avtype == "": for snclq in snclqs: - snclqList = snclq.split('.') - n =snclqList[0] + snclqList = snclq.split(".") + n = snclqList[0] s = snclqList[1] l = snclqList[2] - if l == '': - luse = '--' + if l == "": + luse = "--" else: luse = l c = snclqList[3] q = snclqList[4] - + if q == "M": service = "fdsnws" elif q == "D": service = "ph5ws" - + if service not in services: services.append(service) - URL = f'http://service.iris.edu/{service}/availability/1/query?format=text&' \ - f'net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&' \ - f'starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&' \ - f'mergegaps={tolerance}&includerestricted=true&nodata=404' - try: - tmpDF = pd.read_csv(URL, sep=' ', dtype={'Location': str, 'Station': str}, parse_dates=['Earliest','Latest']) - tmpDF['staloc'] = f'{s}.{luse}' + URL = ( + f"http://service.iris.edu/{service}/availability/1/query?format=text&" + f"net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&" + f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" + f"mergegaps={tolerance}&includerestricted=true&nodata=404" + ) + try: + tmpDF = pd.read_csv( + URL, + sep=" ", + dtype={"Location": str, "Station": str}, + parse_dates=["Earliest", "Latest"], + ) + tmpDF["staloc"] = f"{s}.{luse}" availabilityDF = availabilityDF.append(tmpDF, ignore_index=True) + print(availabilityDF) except Exception as e: pass - - - elif avtype == 'extents': + + elif avtype == "extents": # Don't loop over the stations unless necessary - nets = ','.join(list(set([n.split('.')[0] for n in snclqs]))) - stas = ','.join(list(set([n.split('.')[1] for n in snclqs]))) + nets = ",".join(list(set([n.split(".")[0] for n in snclqs]))) + stas = ",".join(list(set([n.split(".")[1] for n in snclqs]))) stas = "*" - locs = ','.join(list(set([n.split('.')[2] for n in snclqs]))) + locs = ",".join(list(set([n.split(".")[2] for n in snclqs]))) if locs == "": - locs = '*' - chans = ','.join(list(set([n.split('.')[3] for n in snclqs]))) - qs = ','.join(list(set([n.split('.')[4] for n in snclqs]))) - - if qs == 'D': - service = 'ph5ws' - elif qs == 'M': - service = 'fdsnws' - + locs = "*" + chans = ",".join(list(set([n.split(".")[3] for n in snclqs]))) + qs = ",".join(list(set([n.split(".")[4] for n in snclqs]))) + + if qs == "D": + service = "ph5ws" + elif qs == "M": + service = "fdsnws" + if service not in services: services.append(service) - - URL = f'http://service.iris.edu/{service}/availability/1/extent?format=text&' \ - f'net={nets}&sta={stas}&loc={locs}&cha={chans}&quality={qs}&' \ - f'starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&' \ - f'includerestricted=true&nodata=404' + + URL = ( + f"http://service.iris.edu/{service}/availability/1/extent?format=text&" + f"net={nets}&sta={stas}&loc={locs}&cha={chans}&quality={qs}&" + f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" + f"includerestricted=true&nodata=404" + ) try: - availabilityDF = pd.read_csv(URL, sep=' ', dtype={'Location': str, 'Station': str}, parse_dates=['Earliest','Latest']) + availabilityDF = pd.read_csv( + URL, + sep=" ", + dtype={"Location": str, "Station": str}, + parse_dates=["Earliest", "Latest"], + ) except: - print(" INFO: Unable to retrieve availability in one go, trying to loop over stations instead") - + print( + " INFO: Unable to retrieve availability in one go, trying to loop over stations instead" + ) + for snclq in snclqs: - snclqList = snclq.split('.') - n =snclqList[0] + snclqList = snclq.split(".") + n = snclqList[0] s = snclqList[1] l = snclqList[2] - if l == '': - luse = '--' + if l == "": + luse = "--" else: luse = l c = snclqList[3] q = snclqList[4] - + if q == "M": service = "fdsnws" elif q == "D": service = "ph5ws" - + if service not in services: services.append(service) - - URL = f'http://service.iris.edu/{service}/availability/1/extent?format=text&' \ - f'net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&' \ - f'starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&' \ - f'&includerestricted=true&nodata=404' - - try: - tmpDF = pd.read_csv(URL, sep=' ', dtype={'Location': str, 'Station': str}, parse_dates=['Earliest','Latest']) - tmpDF['staloc'] = f'{s}.{luse}' + + URL = ( + f"http://service.iris.edu/{service}/availability/1/extent?format=text&" + f"net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&" + f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" + f"&includerestricted=true&nodata=404" + ) + + try: + tmpDF = pd.read_csv( + URL, + sep=" ", + dtype={"Location": str, "Station": str}, + parse_dates=["Earliest", "Latest"], + ) + tmpDF["staloc"] = f"{s}.{luse}" availabilityDF = availabilityDF.append(tmpDF, ignore_index=True) + print(availabilityDF) except: pass - - - - - - -# availabilityDF = availabilityDF.apply(lambda x: x.str.strip() if x.dtype == "object" else x) + + # availabilityDF = availabilityDF.apply(lambda x: x.str.strip() if x.dtype == "object" else x) availabilityDF.rename(columns=lambda x: x.strip().lower(), inplace=True) - availabilityDF.rename(columns = {'#network': 'network'}, inplace=True) - - return availabilityDF, services - + availabilityDF.rename(columns={"#network": "network"}, inplace=True) + + return availabilityDF, services + + +def getDataStartEnd( + target, + thisNetwork, + thisStation, + thisLocation, + thisChannel, + startDate, + endDate, + tolerance, + availabilityExists, +): + """Get the start/end times of the data for a target, via the availability + service if it's up, otherwise falling back to MUSTANG extents.""" + if availabilityExists: + thisStationAvDF, _ = getAvailability( + [target], startDate, endDate, tolerance, "" + ) + thisAvail = thisStationAvDF[ + (thisStationAvDF["station"] == thisStation) + & (thisStationAvDF["network"] == thisNetwork) + & (thisStationAvDF["location"] == thisLocation) + & (thisStationAvDF["channel"] == thisChannel) + ] + thisStart = thisAvail["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() + thisEnd = thisAvail["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() + else: + thisStationAvDF = get_extents_from_mustang([target], startDate, endDate) + thisStart = thisStationAvDF["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() + thisEnd = thisStationAvDF["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() + + return thisStart, thisEnd + def retrieveMetrics(URL, metric): - response = requests.get(URL) + response = requests.get(URL) tempDF = pd.read_csv(StringIO(response.text), header=1) - tempDF.rename(columns={'target':'snclq'}, inplace=True) - tempDF['target'] = tempDF['snclq'].apply(lambda x: '.'.join(x.split('.')[0:4])) - tempDF['station'] = tempDF['snclq'].apply(lambda x: '.'.join(x.split('.')[1:3])) ## Because "station" is really "station.location" - tempDF['station'] = [x + '--' if x.endswith('.') else x for x in tempDF['station'] ] - - if (not metric == 'transfer_function') and (not metric == 'orientation_check'): - tempDF.rename(columns = {'value': metric}, inplace=True) - tempDF[metric] = tempDF[metric].map(float) - tempDF.drop('lddate', axis=1, inplace=True) - - tempDF['start'] = pd.to_datetime(tempDF['start']) - tempDF['end'] = pd.to_datetime(tempDF['end']) - + tempDF.rename(columns={"target": "snclq"}, inplace=True) + tempDF["target"] = tempDF["snclq"].apply(lambda x: ".".join(x.split(".")[0:4])) + tempDF["station"] = tempDF["snclq"].apply( + lambda x: ".".join(x.split(".")[1:3]) + ) ## Because "station" is really "station.location" + tempDF["station"] = [x + "--" if x.endswith(".") else x for x in tempDF["station"]] + + if (not metric == "transfer_function") and (not metric == "orientation_check"): + tempDF.rename(columns={"value": metric}, inplace=True) + tempDF[metric] = tempDF[metric].map(float) + tempDF.drop("lddate", axis=1, inplace=True) + + tempDF["start"] = pd.to_datetime(tempDF["start"]) + tempDF["end"] = pd.to_datetime(tempDF["end"]) + if metric == "sample_rms" and tempDF.empty: print(f" --> Unable to get measurements for {metric}") print(f" {URL}") - + return tempDF -def addMetricToDF(metric, DF, network, stations, locations, channels, startDate, endDate): - if not (metric== 'ts_percent_availability_total' or metric == 'percent_availability'): + +def checkTsMetricsExist(startDate, endDate, timeout=15): + # Probe against a long-running reference station (IU.ANMO) rather than the + # network being reported on, so the result reflects whether ts_ metrics exist + # at all rather than whether this particular network has data for them. + testURL = ( + f"https://service.earthscope.org/mustang/measurements/1/query?" + f"metric=ts_percent_availability_total&net=IU&sta=ANMO&loc=00&chan=BHZ&" + f"format=text&timewindow={startDate},{endDate}&nodata=404" + ) + try: + response = requests.get(testURL, timeout=timeout) + if "Could not find metric" in response.text: + return False + return True + except Exception: + return False + + +def addMetricToDF( + metric, DF, network, stations, locations, channels, startDate, endDate +): + if not ( + metric == "ts_percent_availability_total" or metric == "percent_availability" + ): print(f" Retrieving {metric}") chanList = list() - for chan in channels.split(','): - if (len(chan) == 3) and (chan[2] in ['*', '?', '.', '_']): - chan = f'{chan[0:2]}Z' + for chan in channels.split(","): + if (len(chan) == 3) and (chan[2] in ["*", "?", ".", "_"]): + chan = f"{chan[0:2]}Z" elif len(chan) == 2: chan = f"{chan}Z" elif chan == "*": chan = "??Z" chanList.append(chan) - URL = f"http://service.iris.edu/mustang/measurements/1/query?metric={metric}&net={network}&" \ - f"sta={stations}&loc={locations}&chan={','.join(chanList)}" \ - f'&format=text&timewindow={startDate},{endDate}&nodata=404' - - + URL = ( + f"http://service.iris.edu/mustang/measurements/1/query?metric={metric}&net={network}&" + f"sta={stations}&loc={locations}&chan={','.join(chanList)}" + f"&format=text&timewindow={startDate},{endDate}&nodata=404" + ) + try: tempDF = retrieveMetrics(URL, metric) except Exception as e: - if not metric== 'ts_percent_availability_total': - print(f" --> Unable to get measurements for {metric}, waiting 5 seconds and trying again") + if not metric == "ts_percent_availability_total": + print( + f" --> Unable to get measurements for {metric}, waiting 5 seconds and trying again" + ) print(f" {URL}") -# print(f" {e}") + # print(f" {e}") time.sleep(5) try: tempDF = retrieveMetrics(URL, metric) except: - if not metric== 'ts_percent_availability_total': - print(f" --> Still unable to get measurements for {metric}, bypassing" ) + if not metric == "ts_percent_availability_total": + print( + f" --> Still unable to get measurements for {metric}, bypassing" + ) tempDF = pd.DataFrame() - - + if DF.empty: DF = tempDF.copy() else: try: - DF = pd.merge(DF, tempDF, how='outer', left_on=['target','snclq', 'station', 'start', 'end'], right_on=['target','snclq','station', 'start', 'end']) + DF = pd.merge( + DF, + tempDF, + how="outer", + left_on=["target", "snclq", "station", "start", "end"], + right_on=["target", "snclq", "station", "start", "end"], + ) except: print(f" ERROR: Something went wrong with the {metric}") - + return DF - + + def getMetadata(network, stations, locations, channels, startDate, endDate, level): # This one and getStations are almost redundant, except that they return at different # levels. Merge into one function that also takes level as an input? - - if level == 'channel': + + if level == "channel": stationDF = pd.DataFrame() - if level == 'station': - stationDF = pd.DataFrame(columns=['#Network', 'Station', 'Latitude' , 'Longitude' , 'Elevation' , 'SiteName' , 'StartTime' , 'EndTime']) - if level == 'network': - stationDF = pd.DataFrame(columns=['#Network',]) - + if level == "station": + stationDF = pd.DataFrame( + columns=[ + "#Network", + "Station", + "Latitude", + "Longitude", + "Elevation", + "SiteName", + "StartTime", + "EndTime", + ] + ) + if level == "network": + stationDF = pd.DataFrame( + columns=[ + "#Network", + ] + ) + chanList = list() - for chan in channels.split(','): - if (len(chan) == 3) and (chan[2] in ['*', '?', '.', '_']): - chan = f'{chan[0:2]}Z' + for chan in channels.split(","): + if (len(chan) == 3) and (chan[2] in ["*", "?", ".", "_"]): + chan = f"{chan[0:2]}Z" elif len(chan) == 2: -# chan = f"{chan}Z,{chan}1" + # chan = f"{chan}Z,{chan}1" chan = f"{chan}Z" elif chan == "*": -# chan = "??Z,??1" + # chan = "??Z,??1" chan = "??Z" chanList.append(chan) - try: - # Call Fed Catalog to know what service the network can be retrieved using. + # Call Fed Catalog to know what service the network can be retrieved using. print(" Calling on Fed Catalog") - fedURL = f"http://service.iris.edu/irisws/fedcatalog/1/query?" \ - f"net={network}&sta={stations}&loc={locations}&cha={','.join(chanList)}&" \ - f"starttime={startDate}&endtime={endDate}" \ - f"&format=request&includeoverlaps=false" - + fedURL = ( + f"http://service.iris.edu/irisws/fedcatalog/1/query?" + f"net={network}&sta={stations}&loc={locations}&cha={','.join(chanList)}&" + f"starttime={startDate}&endtime={endDate}" + f"&format=request&includeoverlaps=false" + ) + try: with urllib.request.urlopen(fedURL) as response: - html_content = response.read().decode('utf-8') - + html_content = response.read().decode("utf-8") + services = [] - for ln in html_content.split('\n'): + for ln in html_content.split("\n"): if ln.startswith("STATIONSERVICE="): - serviceURL = ln.split('=')[1] - if 'iris' in serviceURL: + serviceURL = ln.split("=")[1] + if "earthscope" in serviceURL: services.append(serviceURL) - + except Exception as e: - print(" ERROR: unable to retrieve fed catalog information about where the data lives - %s\n%s " % (fedURL, e)) - services = ['http://service.iris.edu/fdsnws/station/1/', 'http://service.iris.edu/ph5ws/station/1/'] - + print( + " ERROR: unable to retrieve fed catalog information about where the data lives - %s\n%s " + % (fedURL, e) + ) + services = [ + "http://service.iris.edu/fdsnws/station/1/", + "http://service.iris.edu/ph5ws/station/1/", + ] + for service in services: # To prevent needing to know a priori where it's from, try both and only add if attempt is successful # Most experiments are one-archive only, but some have been split in the past try: print(" Calling on Station Service") - stationURL = f"{service}query?" \ - f"net={network}&sta={stations}&loc={locations}&cha={','.join(chanList)}&" \ - f"starttime={startDate}&endtime={endDate}&level={level}" \ - f"&format=text&includecomments=true&nodata=404" + stationURL = ( + f"{service}query?" + f"net={network}&sta={stations}&loc={locations}&cha={','.join(chanList)}&" + f"starttime={startDate}&endtime={endDate}&level={level}" + f"&format=text&includecomments=true&nodata=404" + ) - if level == 'channel': + if level == "channel": try: - tmpDF = pd.read_csv(stationURL, sep='|', dtype={' Location ': str, ' Station ': str}) + tmpDF = pd.read_csv( + stationURL, + sep="|", + dtype={" Location ": str, " Station ": str}, + ) tmpDF.rename(columns=lambda x: x.strip(), inplace=True) - tmpDF.rename(columns = {'#Network': 'Network'}, inplace=True) - tmpDF['Location'] = tmpDF.Location.replace(np.nan, '', regex=True) - tmpDF['Target'] = tmpDF[['Network', 'Station', 'Location','Channel']].apply(lambda x: '.'.join(x.map(str)), axis=1) + tmpDF.rename(columns={"#Network": "Network"}, inplace=True) + tmpDF["Location"] = tmpDF.Location.replace( + np.nan, "", regex=True + ) + tmpDF["Target"] = tmpDF[ + ["Network", "Station", "Location", "Channel"] + ].apply(lambda x: ".".join(x.map(str)), axis=1) tmpDF.columns = tmpDF.columns.str.lower() - - tmpDF['starttime'] = pd.to_datetime(tmpDF['starttime']) - tmpDF['endtime'] = pd.to_datetime(tmpDF['endtime']) - + + tmpDF["starttime"] = pd.to_datetime(tmpDF["starttime"]) + tmpDF["endtime"] = pd.to_datetime(tmpDF["endtime"]) + except Exception as e: - print(f" ERROR: Unable to retrieve channel information from {stationURL}") - - elif level == 'station' or level == 'network': + print( + f" ERROR: Unable to retrieve channel information from {stationURL}" + ) + + elif level == "station" or level == "network": try: - tmpDF = pd.read_csv(stationURL, sep='|') + tmpDF = pd.read_csv(stationURL, sep="|") tmpDF.rename(columns=lambda x: x.strip(), inplace=True) except: - print(f" ERROR: Unable to retrieve metadata information from {stationURL}") + print( + f" ERROR: Unable to retrieve metadata information from {stationURL}" + ) except: tmpDF = pd.DataFrame() - + stationDF = pd.concat([stationDF, tmpDF], ignore_index=True) + print(stationDF) except: - print(" ERROR: Unable to retrieve metadata") + print(" ERROR: Unable to retrieve metadata") return stationDF - + return stationDF + def retrieveExpectedPDFs(NSLC, startDate, endDate): - net = NSLC.split('.')[0] - sta = NSLC.split('.')[1] - loc = NSLC.split('.')[2] + net = NSLC.split(".")[0] + sta = NSLC.split(".")[1] + loc = NSLC.split(".")[2] cha = f"{NSLC.split('.')[3]}?" - URL = f'http://service.iris.edu/mustang/noise-pdf-browser/1/availability?network={net}&station={sta}&location={loc}&channel={cha}&starttime={startDate}&endtime={endDate}&interval=all' + URL = f"http://service.iris.edu/mustang/noise-pdf-browser/1/availability?network={net}&station={sta}&location={loc}&channel={cha}&starttime={startDate}&endtime={endDate}&interval=all" - response = requests.get(URL) + response = requests.get(URL) if response.text.startswith("Error"): # Wait 5 seconds and try again - print(f" --> Error retrieving list of expected PDFs for {NSLC}, waiting 5 seconds and trying again") + print( + f" --> Error retrieving list of expected PDFs for {NSLC}, waiting 5 seconds and trying again" + ) time.sleep(5) - response = requests.get(URL) + response = requests.get(URL) if response.text.startswith("Error"): print(f" --> Unable to retrieve PDF list for {NSLC}") -# print(response.text) + # print(response.text) expectedTargets = list() - - # doing it this way so that this section will run if either the first or second attempt was successful - if not response.text.startswith("Error"): - expectedTargets = [x.split(',')[0] for x in response.text.split('\n') if not x == ''] - + + # doing it this way so that this section will run if either the first or second attempt was successful + if not response.text.startswith("Error"): + expectedTargets = [ + x.split(",")[0] for x in response.text.split("\n") if not x == "" + ] + return expectedTargets + def getPDF(target, startDate, endDate, spectPowerRange, imageDir): - net = target.split('.')[0] - sta = target.split('.')[1] - loc = target.split('.')[2] - if loc == '': - loc = '--' - cha = target.split('.')[3] - - - plot_titlefont=20 - plot_subtitlefont=18 - plot_axisfont=16 - plot_labelfont=18 - plotArguments = f"plot.titlefont.size={plot_titlefont}&plot.subtitlefont.size={plot_subtitlefont}" \ - f"&plot.axisfont.size={plot_axisfont}&plot.labelfont.size={plot_labelfont}" - - URL = f"http://service.iris.edu/mustang/noise-pdf/1/query?&" \ - f"network={net}&station={sta}&location={loc}&channel={cha}&quality=?&" \ - f"starttime={startDate}&endtime={endDate}&format=plot&plot.interpolation=bicubic&nodata=404&" \ - f"plot.power.min={spectPowerRange[0]}&plot.power.max={spectPowerRange[1]}&{plotArguments}" + net = target.split(".")[0] + sta = target.split(".")[1] + loc = target.split(".")[2] + if loc == "": + loc = "--" + cha = target.split(".")[3] + + plot_titlefont = 20 + plot_subtitlefont = 18 + plot_axisfont = 16 + plot_labelfont = 18 + plotArguments = ( + f"plot.titlefont.size={plot_titlefont}&plot.subtitlefont.size={plot_subtitlefont}" + f"&plot.axisfont.size={plot_axisfont}&plot.labelfont.size={plot_labelfont}" + ) + + URL = ( + f"http://service.iris.edu/mustang/noise-pdf/1/query?&" + f"network={net}&station={sta}&location={loc}&channel={cha}&quality=?&" + f"starttime={startDate}&endtime={endDate}&format=plot&plot.interpolation=bicubic&nodata=404&" + f"plot.power.min={spectPowerRange[0]}&plot.power.max={spectPowerRange[1]}&{plotArguments}" + ) response = requests.get(URL) - filename = (f"{imageDir}/{target}_PDF.png").replace('*','').replace('?','') - + filename = (f"{imageDir}/{target}_PDF.png").replace("*", "").replace("?", "") + file = open(filename, "wb") file.write(response.content) file.close() - + return filename - -def getSpectrogram(target, startDate, endDate, spectPowerRange, spectColorPalette, imageDir): - powerRange = ','.join([str(x) for x in spectPowerRange]) - plot_titlefont=20 - plot_subtitlefont=18 - plot_axisfont=16 - plot_labelfont=18 - plotArguments = f"plot.titlefont.size={plot_titlefont}&plot.subtitlefont.size={plot_subtitlefont}" \ - f"&plot.axisfont.size={plot_axisfont}&plot.labelfont.size={plot_labelfont}" - - URL = f"http://service.iris.edu/mustang/noise-spectrogram/1/query?target={target}&" \ - f"starttime={startDate}&endtime={endDate}&output=power&format=plot&plot.color.palette={spectColorPalette}&" \ - f"plot.powerscale.range={powerRange}&plot.horzaxis=time&plot.time.matchrequest=true&{plotArguments}&" \ - f"plot.time.tickunit=auto&plot.time.invert=false&plot.powerscale.show=true&plot.powerscale.orientation=horz&nodata=404" - - + + +def getSpectrogram( + target, startDate, endDate, spectPowerRange, spectColorPalette, imageDir +): + powerRange = ",".join([str(x) for x in spectPowerRange]) + plot_titlefont = 20 + plot_subtitlefont = 18 + plot_axisfont = 16 + plot_labelfont = 18 + plotArguments = ( + f"plot.titlefont.size={plot_titlefont}&plot.subtitlefont.size={plot_subtitlefont}" + f"&plot.axisfont.size={plot_axisfont}&plot.labelfont.size={plot_labelfont}" + ) + + URL = ( + f"http://service.iris.edu/mustang/noise-spectrogram/1/query?target={target}&" + f"starttime={startDate}&endtime={endDate}&output=power&format=plot&plot.color.palette={spectColorPalette}&" + f"plot.powerscale.range={powerRange}&plot.horzaxis=time&plot.time.matchrequest=true&{plotArguments}&" + f"plot.time.tickunit=auto&plot.time.invert=false&plot.powerscale.show=true&plot.powerscale.orientation=horz&nodata=404" + ) + response = requests.get(URL) filename = f"{imageDir}/{target}_spectrogram.png" file = open(filename, "wb") file.write(response.content) file.close() - + return filename def getBoundsZoomLevel(bounds, mapDim): + """ + source: https://stackoverflow.com/questions/6048975/google-maps-v3-how-to-calculate-the-zoom-level-for-a-given-bounds + :param bounds: list of ne and sw lat/lon + :param mapDim: dictionary with image size in pixels + :return: zoom level to fit bounds in the visible area + """ + n_lat = bounds[0] + w_long = bounds[1] + s_lat = bounds[2] + e_long = bounds[3] + + scale = ( + 2 # adjustment to reflect MapBox base tiles are 512x512 vs. Google's 256x256 + ) + WORLD_DIM = {"height": 256 * scale, "width": 256 * scale} + ZOOM_MAX = 16 + ZOOM_MIN = 0.5 + + def latRad(lat): + sin = np.sin(lat * np.pi / 180) + radX2 = np.log((1 + sin) / (1 - sin)) / 2 + return max(min(radX2, np.pi), -np.pi) / 2 + + def zoom(mapPx, worldPx, fraction): + return np.floor(np.log(mapPx / worldPx / fraction) / np.log(2)) + + latFraction = (latRad(n_lat) - latRad(s_lat)) / np.pi + + lngDiff = e_long - w_long + lngFraction = ((lngDiff + 360) if lngDiff < 0 else lngDiff) / 360 + lat_lng_fraction = max(latFraction, lngFraction) / min(latFraction, lngFraction) + + latZoom = zoom(mapDim["height"], WORLD_DIM["height"], latFraction) + lngZoom = zoom(mapDim["width"], WORLD_DIM["width"], lngFraction) + zoomDiff = max(latZoom, lngZoom) + + if lat_lng_fraction > 10: + zoomDiff = max(latZoom, lngZoom) - min(latZoom, lngZoom) + + return min(max(zoomDiff, ZOOM_MIN), ZOOM_MAX) + - """ - source: https://stackoverflow.com/questions/6048975/google-maps-v3-how-to-calculate-the-zoom-level-for-a-given-bounds - :param bounds: list of ne and sw lat/lon - :param mapDim: dictionary with image size in pixels - :return: zoom level to fit bounds in the visible area - """ - n_lat = bounds[0] - w_long = bounds[1] - s_lat = bounds[2] - e_long = bounds[3] - - scale = 2 # adjustment to reflect MapBox base tiles are 512x512 vs. Google's 256x256 - WORLD_DIM = {'height': 256 * scale, 'width': 256 * scale} - ZOOM_MAX = 16 - ZOOM_MIN = 0.5 - - def latRad(lat): - sin = np.sin(lat * np.pi / 180) - radX2 = np.log((1 + sin) / (1 - sin)) / 2 - return max(min(radX2, np.pi), -np.pi) / 2 - - def zoom(mapPx, worldPx, fraction): - return np.floor(np.log(mapPx / worldPx / fraction) / np.log(2)) - - latFraction = (latRad(n_lat) - latRad(s_lat)) / np.pi - - lngDiff = e_long - w_long - lngFraction = ((lngDiff + 360) if lngDiff < 0 else lngDiff) / 360 - lat_lng_fraction = max(latFraction, lngFraction) / min(latFraction, lngFraction) - - latZoom = zoom(mapDim['height'], WORLD_DIM['height'], latFraction) - lngZoom = zoom(mapDim['width'], WORLD_DIM['width'], lngFraction) - zoomDiff = max(latZoom, lngZoom) - - if lat_lng_fraction > 10: - zoomDiff = max(latZoom, lngZoom) - min(latZoom, lngZoom) - - return min(max(zoomDiff, ZOOM_MIN), ZOOM_MAX) - - def getMetricLabel(metric): - metricLabels = {'amplifier_saturation':'daily flag count \n(number of occurrences)', - 'calibration_signal':'daily flag count \n(number of occurrences)', - 'clock_locked':'daily flag count \n(number of occurrences)', - 'cross_talk':'correlation coefficient \n', # no units - 'data_latency':'latency (seconds) \n', - 'dc_offset':'daily indicator of likelihood of \nDC offset shift', # no units - 'dead_channel_gsn':'indicator \n', - 'dead_channel_lin':'standard deviation of residuals (dB) \n', - 'digital_filter_charging':'daily flag count \n(number of occurrences)', - 'digitizer_clipping':'daily flag count \n(number of occurrences)', - 'event_begin':'daily flag count \n(number of occurrences)', - 'event_end':'daily flag count \n(number of occurrences)', - 'event_in_progress':'daily flag count \n(number of occurrences)', - 'feed_latency':'latency (seconds) \n', - 'gap_list':'daily gap length \n(seconds)', - 'glitches':'daily flag count \n(number of occurrences)', - 'max_gap':'daily maximum gap length \n(seconds)', - 'max_overlap':'daily overlap length \n(seconds)', - 'max_range':'daily maximum amplitude range, \nwindowed (counts)', - 'max_stalta':'daily \nshort-term average / long-term \naverage', # no units - 'missing_padded_data':'daily flag count \n(number of occurrences)', - 'num_gaps':'daily gap count \n(number of occurrences)', - 'num_overlaps':'daily overlap count \n(number of occurrences)', - 'num_spikes':'daily outlier count \n(number of occurrences)', - 'pct_above_nhnm':'daily PDF matrix above \nNew High Noise Model (%)', - 'pct_below_nlnm':'daily PDF matrix below \nNew Low Noise Model (%)', - 'percent_availability':'daily availability (%) \n', - 'polarity_check':'maximum cross-correlation \nfunction', # no units - 'pressure_effects':'daily zero-lag \ncross-correlation function', # no units - 'sample_max':'daily maximum amplitude \n(counts)', - 'sample_mean':'daily mean amplitude \n(counts)', - 'sample_median':'daily median amplitude \n(counts)', - 'sample_min':'daily minimum amplitude \n(counts)', - 'sample_rate_channel':'daily indicator \n', - 'sample_rate_resp':'daily indicator \n', - 'sample_rms':'daily root-mean-square variance (counts) \n', - 'scale_corrected_sample_rms':'daily root-mean-squared variance,\nscaled by sensitivity', - 'sample_snr':'signal-to-noise ratio \n', # no units - 'sample_unique':'daily unique sample values \n(number of occurrences)', - 'spikes':'daily flag count \n(number of occurrences)', - 'suspect_time_tag':'daily flag count \n(number of occurrences)', - 'telemetry_sync_error':'daily flag count \n(number of occurrences)', - 'timing_correction':'daily flag count \n(number of occurrences)', - 'timing_quality':'daily average timing quality (%) \n', - 'total_latency':'latency (seconds) \n', - 'ts_num_gaps':'daily gap count \n(number of occurrences)', - 'ts_num_gaps_total':'gap count \n(number of occurrences)', - 'ts_max_gap':'daily maximum gap length \n(seconds)', - 'ts_max_gap_total':'maximum gap length \n(seconds)', - 'ts_gap_length':'daily total gap length \n(seconds)', - 'ts_gap_length_total':'total gap length (seconds) \n', - 'ts_percent_availability':'daily availability (%) \n', - 'ts_percent_availability_total':'availability (%) \n', - 'ts_channel_up_time':'daily trace segment length \n(seconds)', - 'ts_channel_continuity':'trace segment length (seconds) \n', - 'gain_ratio':'data/metadata gain ratio', # no units - 'phase_diff':'data-metadata phase difference \n(degrees)', - 'ms_coherence':'coherence function \n', # no units - } - + metricLabels = { + "amplifier_saturation": "daily flag count \n(number of occurrences)", + "calibration_signal": "daily flag count \n(number of occurrences)", + "clock_locked": "daily flag count \n(number of occurrences)", + "cross_talk": "correlation coefficient \n", # no units + "data_latency": "latency (seconds) \n", + "dc_offset": "daily indicator of likelihood of \nDC offset shift", # no units + "dead_channel_gsn": "indicator \n", + "dead_channel_lin": "standard deviation of residuals (dB) \n", + "digital_filter_charging": "daily flag count \n(number of occurrences)", + "digitizer_clipping": "daily flag count \n(number of occurrences)", + "event_begin": "daily flag count \n(number of occurrences)", + "event_end": "daily flag count \n(number of occurrences)", + "event_in_progress": "daily flag count \n(number of occurrences)", + "feed_latency": "latency (seconds) \n", + "gap_list": "daily gap length \n(seconds)", + "glitches": "daily flag count \n(number of occurrences)", + "max_gap": "daily maximum gap length \n(seconds)", + "max_overlap": "daily overlap length \n(seconds)", + "max_range": "daily maximum amplitude range, \nwindowed (counts)", + "max_stalta": "daily \nshort-term average / long-term \naverage", # no units + "missing_padded_data": "daily flag count \n(number of occurrences)", + "num_gaps": "daily gap count \n(number of occurrences)", + "num_overlaps": "daily overlap count \n(number of occurrences)", + "num_spikes": "daily outlier count \n(number of occurrences)", + "pct_above_nhnm": "daily PDF matrix above \nNew High Noise Model (%)", + "pct_below_nlnm": "daily PDF matrix below \nNew Low Noise Model (%)", + "percent_availability": "daily availability (%) \n", + "polarity_check": "maximum cross-correlation \nfunction", # no units + "pressure_effects": "daily zero-lag \ncross-correlation function", # no units + "sample_max": "daily maximum amplitude \n(counts)", + "sample_mean": "daily mean amplitude \n(counts)", + "sample_median": "daily median amplitude \n(counts)", + "sample_min": "daily minimum amplitude \n(counts)", + "sample_rate_channel": "daily indicator \n", + "sample_rate_resp": "daily indicator \n", + "sample_rms": "daily root-mean-square variance (counts) \n", + "scale_corrected_sample_rms": "daily root-mean-squared variance,\nscaled by sensitivity", + "sample_snr": "signal-to-noise ratio \n", # no units + "sample_unique": "daily unique sample values \n(number of occurrences)", + "spikes": "daily flag count \n(number of occurrences)", + "suspect_time_tag": "daily flag count \n(number of occurrences)", + "telemetry_sync_error": "daily flag count \n(number of occurrences)", + "timing_correction": "daily flag count \n(number of occurrences)", + "timing_quality": "daily average timing quality (%) \n", + "total_latency": "latency (seconds) \n", + "ts_num_gaps": "daily gap count \n(number of occurrences)", + "ts_num_gaps_total": "gap count \n(number of occurrences)", + "ts_max_gap": "daily maximum gap length \n(seconds)", + "ts_max_gap_total": "maximum gap length \n(seconds)", + "ts_gap_length": "daily total gap length \n(seconds)", + "ts_gap_length_total": "total gap length (seconds) \n", + "ts_percent_availability": "daily availability (%) \n", + "ts_percent_availability_total": "availability (%) \n", + "ts_channel_up_time": "daily trace segment length \n(seconds)", + "ts_channel_continuity": "trace segment length (seconds) \n", + "gain_ratio": "data/metadata gain ratio", # no units + "phase_diff": "data-metadata phase difference \n(degrees)", + "ms_coherence": "coherence function \n", # no units + } + labelText = metricLabels[metric] return labelText From 62cf4d9d06d840d8479f500a76c81e481853422e Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Mon, 27 Jul 2026 15:53:55 -0700 Subject: [PATCH 02/21] if availability and percent_avail are missing, use report dates --- PIQQA.py | 8 ++++---- reportUtils.py | 20 +++++++++++++++++++- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 6e06f10..f69cc7e 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -1240,6 +1240,7 @@ def doPDFs( availabilityExists, ) + foundit = False print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" ) @@ -1248,7 +1249,6 @@ def doPDFs( thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd ) - foundit = False if thisPctAvail > 75: # Then we can use it and break the cycle. foundit = True @@ -1313,6 +1313,7 @@ def doPDFs( availabilityExists, ) + foundit = False print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" ) @@ -1321,7 +1322,6 @@ def doPDFs( thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd ) - foundit = False if thisPctAvail > 75: # Needs to have >75% and not be the same target that we've already selected for lowest RMS # Then we can use it and break the cycle. Otherwise break with foundit still = False so it will use the original choice @@ -1489,6 +1489,7 @@ def doSpectrograms( availabilityExists, ) + foundit = False print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" ) @@ -1497,7 +1498,6 @@ def doSpectrograms( thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd ) - foundit = False if thisPctAvail > 75: # Then we can use it and break the cycle. foundit = True @@ -1561,6 +1561,7 @@ def doSpectrograms( availabilityExists, ) + foundit = False print( f" INFO: Checking percent availability for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} to see if there is as least 75% availability" ) @@ -1569,7 +1570,6 @@ def doSpectrograms( thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd ) - foundit = False if thisPctAvail > 75: # Needs to have >75% and not be the same target that we've already selected for lowest RMS # Then we can use it and break the cycle. Otherwise break with foundit still = False so it will use the original choice diff --git a/reportUtils.py b/reportUtils.py index 0b8449a..597ab9c 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -217,21 +217,39 @@ def getDataStartEnd( availabilityExists, ): """Get the start/end times of the data for a target, via the availability - service if it's up, otherwise falling back to MUSTANG extents.""" + service if it's up, otherwise falling back to MUSTANG extents. Falls back + to the report's own startDate/endDate if no start/end can be determined + either way (e.g. percent_availability is also down or has no data).""" + reportStart = f"{startDate}T00:00:00" + reportEnd = f"{endDate}T00:00:00" + if availabilityExists: thisStationAvDF, _ = getAvailability( [target], startDate, endDate, tolerance, "" ) + if thisStationAvDF.empty: + return reportStart, reportEnd + thisAvail = thisStationAvDF[ (thisStationAvDF["station"] == thisStation) & (thisStationAvDF["network"] == thisNetwork) & (thisStationAvDF["location"] == thisLocation) & (thisStationAvDF["channel"] == thisChannel) ] + if thisAvail.empty: + return reportStart, reportEnd + thisStart = thisAvail["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() thisEnd = thisAvail["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() else: thisStationAvDF = get_extents_from_mustang([target], startDate, endDate) + if thisStationAvDF.empty: + print( + f" WARNING: Unable to determine data start/end for {target} " + f"(percent_availability returned no data from Mustang), using report timespan instead" + ) + return reportStart, reportEnd + thisStart = thisStationAvDF["earliest"].dt.strftime("%Y-%m-%dT%H:%M:%S").min() thisEnd = thisStationAvDF["latest"].dt.strftime("%Y-%m-%dT%H:%M:%S").max() From 009e97e223143e0ce30eee54cd3269219ff2c342 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Mon, 27 Jul 2026 16:36:30 -0700 Subject: [PATCH 03/21] fix missing spectrograms --- PIQQA.py | 46 +++++++++++++++++++++++++++++----------------- 1 file changed, 29 insertions(+), 17 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index f69cc7e..56b0977 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -128,27 +128,27 @@ def doAvailability( imageDir, tsMetricsExist, ): - # test out whether the availability service is up and running, if not, then skip the availability plot + # test out whether the availability service is up and running; if not, fall back to + # MUSTANG percent_availability extents (no gap-level detail) rather than skipping the plot availability_exists = True availablity_test_url = "https://service.earthscope.org/fdsnws/availability/1" - # if the service is returning an HTTP 410, then bypass the availablity plot try: urllib.request.urlopen(availablity_test_url) except urllib.error.HTTPError as e: if e.code == 410: - print("WARNING: Availability service is down, skipping availability plot") + print( + "WARNING: Availability service is down, falling back to MUSTANG percent_availability for data extents" + ) else: print( - f"WARNING: Availability service returned HTTP {e.code}, skipping availability plot" + f"WARNING: Availability service returned HTTP {e.code}, falling back to MUSTANG percent_availability for data extents" ) availability_exists = False - return {}, splitPlots, {}, [], "", availability_exists except Exception as e: print( - f"WARNING: Availability service is down ({e}), skipping availability plot" + f"WARNING: Availability service is down ({e}), falling back to MUSTANG percent_availability for data extents" ) availability_exists = False - return {}, splitPlots, {}, [], "", availability_exists print("INFO: Producing availability plot") avFilesDictionary = {} # collects plot filenames @@ -290,6 +290,20 @@ def doAvailability( topStations = pctAvDF.iloc[:nTop, :] bottomStations = pctAvDF.iloc[-nBottom:, :] + # Station selection only needs MUSTANG percent-availability data, so this stays + # populated (for PDFs/Spectrograms downstream) even if the availability service is down. + try: + topStationList = topStations["station"].tolist() + except: + topStationList = tmpMetadataDF["station"].tolist()[:nTop] + + topStationsDict[channelGroup] = topStationList + + if not availability_exists: + # Can't draw the extents/gaps plot without the availability service - skip the plot, + # but topStationsDict above still lets PDFs/Spectrograms pick stations for this channel group. + continue + height = max(min(0.3 * nsta, 0.3 * nBoxPlotSta), 2) width = 15 f, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(width, height)) @@ -309,18 +323,11 @@ def doAvailability( if service not in services_used: services_used.append(service) - try: - topStationList = topStations["station"].tolist() - except: - topStationList = tmpMetadataDF["station"].tolist()[:nTop] - topLabels = [ f"{a} ({b:.3f}%)" for a, b in zip(topStations["station"], topStations["availability"]) ] - topStationsDict[channelGroup] = topStationList - datalines = [] metadatalines = [] gaplines = [] @@ -541,7 +548,7 @@ def doAvailability( ax2.add_collection(GapLines) ax2.autoscale() ax2.invert_yaxis() - ax2.set_yticks(np.arange(nTop)) + ax2.set_yticks(np.arange(nBottom)) ax2.set_yticklabels(bottomLabels) ax2.set_xlim( [mpl.dates.datestr2num(startDate), mpl.dates.datestr2num(endDate)] @@ -578,6 +585,13 @@ def doAvailability( ## Then repeat for the case where they aren't broken up by top/bottom else: + allStationList = pctAvDF["station"].tolist() + topStationsDict[channelGroup] = allStationList + + if not availability_exists: + # Can't draw the extents/gaps plot without the availability service - skip the plot, + # but topStationsDict above still lets PDFs/Spectrograms pick stations for this channel group. + continue height = max(min(0.3 * nsta, 0.3 * nBoxPlotSta), 2) width = 15 @@ -590,12 +604,10 @@ def doAvailability( if service not in services_used: services_used.append(service) - allStationList = pctAvDF["station"].tolist() allLabels = [ f"{a} ({b:.3f}%)" for a, b in zip(pctAvDF["station"], pctAvDF["availability"]) ] - topStationsDict[channelGroup] = allStationList datalines = [] metadatalines = [] From 70a0aba77c0dc6b3147933b14c298d3b05b3b7ca Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Mon, 27 Jul 2026 16:46:19 -0700 Subject: [PATCH 04/21] add user agent to service calls --- reportUtils.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/reportUtils.py b/reportUtils.py index 597ab9c..027a3de 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -27,7 +27,19 @@ import pandas as pd import numpy as np import time -import urllib +import urllib.request + +USER_AGENT = "earthscope-piqqa" + +# pandas' read_csv falls back to urllib for plain http(s) URLs (no fsspec installed), +# so installing a default opener here also sets the User-Agent for every pd.read_csv(url) +# call and any bare urllib.request.urlopen() call in this process, not just SESSION.get(). +_opener = urllib.request.build_opener() +_opener.addheaders = [("User-Agent", USER_AGENT)] +urllib.request.install_opener(_opener) + +SESSION = requests.Session() +SESSION.headers.update({"User-Agent": USER_AGENT}) def get_extents_from_mustang(snclqs, startDate, endDate): @@ -257,7 +269,7 @@ def getDataStartEnd( def retrieveMetrics(URL, metric): - response = requests.get(URL) + response = SESSION.get(URL) tempDF = pd.read_csv(StringIO(response.text), header=1) tempDF.rename(columns={"target": "snclq"}, inplace=True) tempDF["target"] = tempDF["snclq"].apply(lambda x: ".".join(x.split(".")[0:4])) @@ -291,7 +303,7 @@ def checkTsMetricsExist(startDate, endDate, timeout=15): f"format=text&timewindow={startDate},{endDate}&nodata=404" ) try: - response = requests.get(testURL, timeout=timeout) + response = SESSION.get(testURL, timeout=timeout) if "Could not find metric" in response.text: return False return True @@ -492,14 +504,14 @@ def retrieveExpectedPDFs(NSLC, startDate, endDate): cha = f"{NSLC.split('.')[3]}?" URL = f"http://service.iris.edu/mustang/noise-pdf-browser/1/availability?network={net}&station={sta}&location={loc}&channel={cha}&starttime={startDate}&endtime={endDate}&interval=all" - response = requests.get(URL) + response = SESSION.get(URL) if response.text.startswith("Error"): # Wait 5 seconds and try again print( f" --> Error retrieving list of expected PDFs for {NSLC}, waiting 5 seconds and trying again" ) time.sleep(5) - response = requests.get(URL) + response = SESSION.get(URL) if response.text.startswith("Error"): print(f" --> Unable to retrieve PDF list for {NSLC}") # print(response.text) @@ -538,7 +550,7 @@ def getPDF(target, startDate, endDate, spectPowerRange, imageDir): f"plot.power.min={spectPowerRange[0]}&plot.power.max={spectPowerRange[1]}&{plotArguments}" ) - response = requests.get(URL) + response = SESSION.get(URL) filename = (f"{imageDir}/{target}_PDF.png").replace("*", "").replace("?", "") file = open(filename, "wb") @@ -568,7 +580,7 @@ def getSpectrogram( f"plot.time.tickunit=auto&plot.time.invert=false&plot.powerscale.show=true&plot.powerscale.orientation=horz&nodata=404" ) - response = requests.get(URL) + response = SESSION.get(URL) filename = f"{imageDir}/{target}_spectrogram.png" file = open(filename, "wb") file.write(response.content) From f53aa7185cacf77d7a4e2c6457670f8b8f7f949a Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Mon, 27 Jul 2026 16:49:11 -0700 Subject: [PATCH 05/21] remove extra print --- reportUtils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/reportUtils.py b/reportUtils.py index 027a3de..8f82df5 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -489,7 +489,6 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve tmpDF = pd.DataFrame() stationDF = pd.concat([stationDF, tmpDF], ignore_index=True) - print(stationDF) except: print(" ERROR: Unable to retrieve metadata") return stationDF From 10adfa073a0468151132de0a2ef22f52b9fca868 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 07:46:38 -0700 Subject: [PATCH 06/21] update service url --> earthscope.org --- PIQQA.py | 38 +++++++++++++++++++------------------- reportUtils.py | 20 ++++++++++---------- 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 56b0977..f731c60 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -1772,7 +1772,7 @@ def doReport( f"This report is intended as a quick, broad overview of the quality of the data archived for the network " f"and time period specified above. IRIS' goal in generating these reports is to give PIs for temporary " f"experiments better insight into the quality of their data archived at the DMC, as well as to demonstrate " - f'the utility provided by the IRIS DMC quality assurance system MUSTANG ' + f'the utility provided by the IRIS DMC quality assurance system MUSTANG ' f"and the many metrics and products it generates. " f"For PIs, we hope that these reports will " f"inform the range of data availability, continuity, and noise levels existing across the network and potentially " @@ -1794,8 +1794,8 @@ def doReport( ###### AVAILABILITY ###### # for service in services: - # stationServiceLink = f'{service}-station service,' - # availabilityServiceLink = f'{service}-availability service,' + # stationServiceLink = f'{service}-station service,' + # availabilityServiceLink = f'{service}-availability service,' # stationServiceLink = f"{stationServiceLink[:-1]}" if availabilityExists: availabilityServiceLink = ( @@ -1818,7 +1818,7 @@ def doReport( f"are from the IRIS DMC {availabilityServiceLink}. " f"The gap tolerance used for this report is {tolerance} seconds.

" f"In addition, each channel lists the percent of data available, using the " - f'{availabilityType} metric' ) # the punctuation at the end of this sentence depends on the metric @@ -1855,7 +1855,7 @@ def doReport( for service in services: serviceLink = ( - f"http://service.iris.edu/{service}/availability/1/query?format=text&" + f"http://service.earthscope.org/{service}/availability/1/query?format=text&" f'net={network}&sta={stations}&loc={locations}&cha={",".join(chans)}&' f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" f"mergegaps={tolerance}&includerestricted=true&nodata=404" @@ -1867,7 +1867,7 @@ def doReport( for service in services: serviceLink = ( - f"http://service.iris.edu/{service}/station/1/query?" + f"http://service.earthscope.org/{service}/station/1/query?" f'net={network}&sta={stations}&loc={locations}&cha={",".join(chans)}&' f"starttime={startDate}&endtime={endDate}&level=channel&" f"format=text&includecomments=true&nodata=404" @@ -1882,7 +1882,7 @@ def doReport( defaultMetricText = "The metrics presented here are" for metric in metricList: defaultMetricText = ( - f'{defaultMetricText} {metric},' ) if metric == metricList[-2]: @@ -1975,15 +1975,15 @@ def doReport( for metric in metricsWithPlots: if metric == "scale_corrected_sample_rms": metric = "sample_rms" - boxPlotIntroText3 = f'{boxPlotIntroText3}
  • {metric}
    • \n' + boxPlotIntroText3 = f'{boxPlotIntroText3}
  • {metric}
    • \n' boxPlotIntroText3 = f"{boxPlotIntroText3} " ###### PDFS ###### pdfPlotIntroText = ( f"

    This section contains Probability Density Function (PDF) plots calculated from Power Spectral Densities (PSD). The first two plots display PDFs " - f'for the stations with the highest and lowest median daily-RMS variance (sample_rms, scaled by metadata sensitivity), ' - f'with an additional criteria that omits station-channels with less than 75% data availability (sample_rms, scaled by metadata sensitivity), ' + f'with an additional criteria that omits station-channels with less than 75% data availability ({availabilityType} metric, ' f"calculated between the first and last timestamp of the data for that channel), unless none of the stations have at least 75% availability. " f"These plots are intended to illustrate the range of characteristic noise levels across the experiment, avoiding stations with significant " @@ -1991,7 +1991,7 @@ def doReport( f"errors. The third PDF is a composite plot of all stations for each channel set, giving an overview of the most common noise levels for the " f"experiment as a whole.

    " f"

    Detailed information about these PDF plots and how MUSTANG generates them " - f'can be found by visiting the ' + f'can be found by visiting the ' f"noise-pdf web service.

    " f"Please note that you can also click the PDF Browser links below to access the MUSTANG PDF Browser, which allows you to view PDF plots " f"for all stations in the network on a total, annual, monthly, daily, or custom time period basis. Within the Browser, you can " @@ -2002,7 +2002,7 @@ def doReport( spectrgramPlotIntoText = ( f"

    The daily-PDF-mode spectrogram plots in this section illustrate the power spectra values across time for two stations with high data " f"availability and continuity. For the purpose of showing potentially different station behavior, they are secondarily selected for " - f'differing median daily-RMS variance (sample_rms ' + f'differing median daily-RMS variance (sample_rms ' f" metric, scaled by metadata sensitivity). For each channel group, the station selection steps are: " f"

      " ) @@ -2019,7 +2019,7 @@ def doReport( if splitPlots == 1: spectrgramPlotIntoText = ( f"{spectrgramPlotIntoText}
    1. Take the list of the {nTop} stations with the highest percent of data availability calculated over the entire requested " - f'report timespan (using {availabilityType}, {spectAvailDescription}).
      2. ' ) else: @@ -2044,7 +2044,7 @@ def doReport( spectrgramPlotIntoText = ( f"{spectrgramPlotIntoText}

      Detailed information about these spectrogram plots and how MUSTANG " - f'generates them can be found by visiting the noise-spectrogram web service.

      ' f"

      Please note that you can also click the Spectrogram Browser links below to access the MUSTANG spectrogram " f"browser for each channel group below, which allows you to view spectrogram plots from all the stations in your network.

      " @@ -2324,7 +2324,7 @@ def doReport( "

      MUSTANG is the Quality Assurance system at the IRIS DMC. It contains around 45 metrics related to the quality of data in the archives there.\n\n" ) f.write( - 'The majority of metrics are available via the measurements web service.\n\n' + 'The majority of metrics are available via the measurements web service.\n\n' ) f.write( 'To learn more about the metrics, navigate to the measurements service Service Interface page and hit the red "Current List of all metrics" for a brief description and links to more detailed documentation.\n\n' @@ -2343,7 +2343,7 @@ def doReport( if metric == "scale_corrected_sample_rms": metric = "sample_rms" metricsLinks = ( - f"http://service.iris.edu/mustang/measurements/1/query?metric={metric}" + f"http://service.earthscope.org/mustang/measurements/1/query?metric={metric}" f"&network={network}&station={stations}&location={locations}&channel={channels}" f"&start={startDate}&end={endDate}&format=text" ) @@ -2381,7 +2381,7 @@ def doReport( nStations = 1 pdfLink = ( - f"http://service.iris.edu/mustang/noise-pdf-browser/1/gallery?" + f"http://service.earthscope.org/mustang/noise-pdf-browser/1/gallery?" f"network={network}&channel={channel[0:2]}?&interval=all&" f"starttime={startDate}&endtime={endDate}" ) @@ -2468,7 +2468,7 @@ def doReport( spectPowerRange = ",".join([str(int) for int in userDefinedPowerRange]) spectLink = ( - f"http://service.iris.edu/mustang/noise-pdf-browser/1/spectrogram?" + f"http://service.earthscope.org/mustang/noise-pdf-browser/1/spectrogram?" f"network={network}&channel={channel[0:2]}?&" f"starttime={startDate}&endtime={endDate}&color.palette={spectColorPalette}&powerrange={spectPowerRange}" ) @@ -2644,7 +2644,7 @@ def main(): --metric[s]=: comma-separated list of metrics to run for the boxplots; defaults: {','.join(metricList)} --maxplot=: number of stations to include in the boxplots; defaults to {nBoxPlotSta} --colorpalette=: color palette for spectrograms; defaults to '{spectColorPalette}' - options available at http://service.iris.edu/mustang/noise-spectrogram/1/ + options available at http://service.earthscope.org/mustang/noise-spectrogram/1/ --includeoutliers=: whether to include outliers in the boxplots, True/False; defaults to {includeOutliers} --spectralrange=: power range to use in the PDFs and spectrograms, comma separated values: min,max; defaults depend on channel type --basemap=: the name of the basemap to be used for the map; defaults to '{basemap}' diff --git a/reportUtils.py b/reportUtils.py index 8f82df5..79bc5af 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -115,7 +115,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): services.append(service) URL = ( - f"http://service.iris.edu/{service}/availability/1/query?format=text&" + f"http://service.earthscope.org/{service}/availability/1/query?format=text&" f"net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&" f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" f"mergegaps={tolerance}&includerestricted=true&nodata=404" @@ -153,7 +153,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): services.append(service) URL = ( - f"http://service.iris.edu/{service}/availability/1/extent?format=text&" + f"http://service.earthscope.org/{service}/availability/1/extent?format=text&" f"net={nets}&sta={stas}&loc={locs}&cha={chans}&quality={qs}&" f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" f"includerestricted=true&nodata=404" @@ -191,7 +191,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): services.append(service) URL = ( - f"http://service.iris.edu/{service}/availability/1/extent?format=text&" + f"http://service.earthscope.org/{service}/availability/1/extent?format=text&" f"net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&" f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" f"&includerestricted=true&nodata=404" @@ -329,7 +329,7 @@ def addMetricToDF( chanList.append(chan) URL = ( - f"http://service.iris.edu/mustang/measurements/1/query?metric={metric}&net={network}&" + f"http://service.earthscope.org/mustang/measurements/1/query?metric={metric}&net={network}&" f"sta={stations}&loc={locations}&chan={','.join(chanList)}" f"&format=text&timewindow={startDate},{endDate}&nodata=404" ) @@ -412,7 +412,7 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve # Call Fed Catalog to know what service the network can be retrieved using. print(" Calling on Fed Catalog") fedURL = ( - f"http://service.iris.edu/irisws/fedcatalog/1/query?" + f"http://service.earthscope.org/irisws/fedcatalog/1/query?" f"net={network}&sta={stations}&loc={locations}&cha={','.join(chanList)}&" f"starttime={startDate}&endtime={endDate}" f"&format=request&includeoverlaps=false" @@ -435,8 +435,8 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve % (fedURL, e) ) services = [ - "http://service.iris.edu/fdsnws/station/1/", - "http://service.iris.edu/ph5ws/station/1/", + "http://service.earthscope.org/fdsnws/station/1/", + "http://service.earthscope.org/ph5ws/station/1/", ] for service in services: @@ -501,7 +501,7 @@ def retrieveExpectedPDFs(NSLC, startDate, endDate): sta = NSLC.split(".")[1] loc = NSLC.split(".")[2] cha = f"{NSLC.split('.')[3]}?" - URL = f"http://service.iris.edu/mustang/noise-pdf-browser/1/availability?network={net}&station={sta}&location={loc}&channel={cha}&starttime={startDate}&endtime={endDate}&interval=all" + URL = f"http://service.earthscope.org/mustang/noise-pdf-browser/1/availability?network={net}&station={sta}&location={loc}&channel={cha}&starttime={startDate}&endtime={endDate}&interval=all" response = SESSION.get(URL) if response.text.startswith("Error"): @@ -543,7 +543,7 @@ def getPDF(target, startDate, endDate, spectPowerRange, imageDir): ) URL = ( - f"http://service.iris.edu/mustang/noise-pdf/1/query?&" + f"http://service.earthscope.org/mustang/noise-pdf/1/query?&" f"network={net}&station={sta}&location={loc}&channel={cha}&quality=?&" f"starttime={startDate}&endtime={endDate}&format=plot&plot.interpolation=bicubic&nodata=404&" f"plot.power.min={spectPowerRange[0]}&plot.power.max={spectPowerRange[1]}&{plotArguments}" @@ -573,7 +573,7 @@ def getSpectrogram( ) URL = ( - f"http://service.iris.edu/mustang/noise-spectrogram/1/query?target={target}&" + f"http://service.earthscope.org/mustang/noise-spectrogram/1/query?target={target}&" f"starttime={startDate}&endtime={endDate}&output=power&format=plot&plot.color.palette={spectColorPalette}&" f"plot.powerscale.range={powerRange}&plot.horzaxis=time&plot.time.matchrequest=true&{plotArguments}&" f"plot.time.tickunit=auto&plot.time.invert=false&plot.powerscale.show=true&plot.powerscale.orientation=horz&nodata=404" From d2fa622efff415c41aeadfc25226dd196643cecb Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 08:06:57 -0700 Subject: [PATCH 07/21] add / to availability url --- PIQQA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PIQQA.py b/PIQQA.py index f731c60..d1aa95e 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -131,7 +131,7 @@ def doAvailability( # test out whether the availability service is up and running; if not, fall back to # MUSTANG percent_availability extents (no gap-level detail) rather than skipping the plot availability_exists = True - availablity_test_url = "https://service.earthscope.org/fdsnws/availability/1" + availablity_test_url = "https://service.earthscope.org/fdsnws/availability/1/" try: urllib.request.urlopen(availablity_test_url) except urllib.error.HTTPError as e: From d68d0a95f49f0dcfeceb00744a9b0593a32fd262 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 08:08:13 -0700 Subject: [PATCH 08/21] update url in readme --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 9f91b28..f918e3d 100644 --- a/README.md +++ b/README.md @@ -122,7 +122,7 @@ Running PIQQA.py without any arguments, or with the -h/--h flags, will produce s --metrics=: comma-separated list of metrics to run for the boxplots; defaults: ts_channel_continuity,sample_rms,ts_num_gaps --maxplot=: number of stations to include in the boxplots; defaults to 30 --colorpalette=: color palette for spectrograms; defaults to 'RdYlBu' - options available at http://service.iris.edu/mustang/noise-spectrogram/1/ + options available at http://service.earthscope.org/mustang/noise-spectrogram/1/ --includeoutliers=: whether to include outliers in the boxplots, True/False; defaults to False --spectralrange=: power range to use in the PDFs and spectrograms, comma separated values: min, max; defaults depend on channel type --basemap=: the name of the basemap to be used for the map; defaults to 'stamen-terrain' @@ -149,7 +149,7 @@ In addition to those fields, there are a number of optional fields. `maxplot` is how you can limit the number of stations displayed in the boxplots. -`colorpalette` determines the color palette used in the PDFs and spectrograms, and a list of all options can be found in the service interface page for the [noise-spectrgram](http://service.iris.edu/mustang/noise-spectrogram/1/) service +`colorpalette` determines the color palette used in the PDFs and spectrograms, and a list of all options can be found in the service interface page for the [noise-spectrgram](http://service.earthscope.org/mustang/noise-spectrogram/1/) service `includeoutliers` toggles the outliers on/off in the boxplots. Generally, it is easier to view the boxplots with the outliers turned off, but there may be cases where they should be displayed. From bc3b9989af26f71b55fe1105a6c70b4631a4c600 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 09:11:02 -0700 Subject: [PATCH 09/21] fix warning message --- PIQQA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PIQQA.py b/PIQQA.py index d1aa95e..6d05643 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -104,7 +104,7 @@ def checkAvailability( except: thisPctAvail = 0 print( - f"**** WARNING: No Percent Availability found for {thisNetwork.thisStation.thisLocation.thisChannel} - are services down? ****" + f"**** WARNING: No Percent Availability found for {thisNetwork}.{thisStation}.{thisLocation}.{thisChannel} - are services down? ****" ) return thisPctAvail From 96d8eda1cc3d0579080cad10652c1b8a56dc2bfb Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 09:31:30 -0700 Subject: [PATCH 10/21] fix links in report --- PIQQA.py | 60 +++++++++++++++++++++++++++++++++----------------------- 1 file changed, 36 insertions(+), 24 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 6d05643..66f4727 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -1798,9 +1798,7 @@ def doReport( # availabilityServiceLink = f'{service}-availability service,' # stationServiceLink = f"{stationServiceLink[:-1]}" if availabilityExists: - availabilityServiceLink = ( - f"http://service.earthscope.org/fdsnws/availability/1/" - ) + availabilityServiceLink = f'fdsnws-availability service' else: availabilityServiceLink = f'mustang service' stationServiceLink = f'fdsnws-station service' @@ -1850,33 +1848,47 @@ def doReport( f"{availabilityIntroText}
      Displaying all stations in the network.
      " ) - try: - availabilityOutroText = "To view the availability numbers used to create the availability plot(s), see:" + if not availabilityExists: + availabilityIntroText = ( + f"{availabilityIntroText}

      " + f"Note: the fdsnws-availability service was unavailable when this report was generated, " + f"so the data-extent and gap plot(s) below could not be produced. Station selection above was " + f"still based on percent-availability metrics from MUSTANG.
      " + ) - for service in services: - serviceLink = ( - f"http://service.earthscope.org/{service}/availability/1/query?format=text&" - f'net={network}&sta={stations}&loc={locations}&cha={",".join(chans)}&' - f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" - f"mergegaps={tolerance}&includerestricted=true&nodata=404" - ) + if services: + try: + availabilityOutroText = "To view the availability numbers used to create the availability plot(s), see:" + + for service in services: + serviceLink = ( + f"http://service.earthscope.org/{service}/availability/1/query?format=text&" + f'net={network}&sta={stations}&loc={locations}&cha={",".join(chans)}&' + f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" + f"mergegaps={tolerance}&includerestricted=true&nodata=404" + ) - availabilityOutroText = f'
      {availabilityOutroText}
      {serviceLink}' + availabilityOutroText = f'
      {availabilityOutroText}
      {serviceLink}' - availabilityOutroText = f"{availabilityOutroText}

      To view the channel metadata time extents used, see: " + availabilityOutroText = f"{availabilityOutroText}

      To view the channel metadata time extents used, see: " - for service in services: - serviceLink = ( - f"http://service.earthscope.org/{service}/station/1/query?" - f'net={network}&sta={stations}&loc={locations}&cha={",".join(chans)}&' - f"starttime={startDate}&endtime={endDate}&level=channel&" - f"format=text&includecomments=true&nodata=404" - ) + for service in services: + serviceLink = ( + f"http://service.earthscope.org/{service}/station/1/query?" + f'net={network}&sta={stations}&loc={locations}&cha={",".join(chans)}&' + f"starttime={startDate}&endtime={endDate}&level=channel&" + f"format=text&includecomments=true&nodata=404" + ) - availabilityOutroText = f'{availabilityOutroText}
      {serviceLink}' + availabilityOutroText = f'{availabilityOutroText}
      {serviceLink}' - except: - availabilityOutroText = "No availability was found for the selected channels." + except: + availabilityOutroText = "No availability was found for the selected channels." + else: + availabilityOutroText = ( + "The availability service was unavailable when this report was generated, " + "so there are no availability/metadata query links to display for this section." + ) ###### BOXPLOTS ###### defaultMetricText = "The metrics presented here are" From 94632076138285694c12db5aefff8d3fa35c289d Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 09:31:48 -0700 Subject: [PATCH 11/21] modernize append->concat for dfs --- reportUtils.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/reportUtils.py b/reportUtils.py index 79bc5af..70e3760 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -128,8 +128,9 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): parse_dates=["Earliest", "Latest"], ) tmpDF["staloc"] = f"{s}.{luse}" - availabilityDF = availabilityDF.append(tmpDF, ignore_index=True) - print(availabilityDF) + availabilityDF = pd.concat( + [availabilityDF, tmpDF], ignore_index=True + ) except Exception as e: pass @@ -205,8 +206,9 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): parse_dates=["Earliest", "Latest"], ) tmpDF["staloc"] = f"{s}.{luse}" - availabilityDF = availabilityDF.append(tmpDF, ignore_index=True) - print(availabilityDF) + availabilityDF = pd.concat( + [availabilityDF, tmpDF], ignore_index=True + ) except: pass From 46551cff235cbf8f095b802ab8c77913ef8c5ec5 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 09:38:55 -0700 Subject: [PATCH 12/21] bump version --- PIQQA.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 66f4727..35d604f 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -24,7 +24,7 @@ """ -version = "v1.0.2" +version = "v1.1.0" import reportUtils import pandas as pd @@ -549,7 +549,7 @@ def doAvailability( ax2.autoscale() ax2.invert_yaxis() ax2.set_yticks(np.arange(nBottom)) - ax2.set_yticklabels(bottomLabels) + ax2.set_yticklabels(bottfomLabels) ax2.set_xlim( [mpl.dates.datestr2num(startDate), mpl.dates.datestr2num(endDate)] ) @@ -1883,7 +1883,9 @@ def doReport( availabilityOutroText = f'{availabilityOutroText}
      {serviceLink}' except: - availabilityOutroText = "No availability was found for the selected channels." + availabilityOutroText = ( + "No availability was found for the selected channels." + ) else: availabilityOutroText = ( "The availability service was unavailable when this report was generated, " From 4367ceffa11b19133dbe82bcdee9e6e7920ae516 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 11:22:24 -0700 Subject: [PATCH 13/21] add handling for when station service fails --- reportUtils.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/reportUtils.py b/reportUtils.py index 70e3760..9a3b49a 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -444,6 +444,7 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve for service in services: # To prevent needing to know a priori where it's from, try both and only add if attempt is successful # Most experiments are one-archive only, but some have been split in the past + tmpDF = pd.DataFrame() try: print(" Calling on Station Service") stationURL = ( @@ -477,6 +478,7 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve print( f" ERROR: Unable to retrieve channel information from {stationURL}" ) + tmpDF = pd.DataFrame() elif level == "station" or level == "network": try: @@ -486,11 +488,17 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve print( f" ERROR: Unable to retrieve metadata information from {stationURL}" ) + tmpDF = pd.DataFrame() except: tmpDF = pd.DataFrame() - stationDF = pd.concat([stationDF, tmpDF], ignore_index=True) + if tmpDF.empty: + continue + if stationDF.empty: + stationDF = tmpDF.copy() + else: + stationDF = pd.concat([stationDF, tmpDF], ignore_index=True) except: print(" ERROR: Unable to retrieve metadata") return stationDF From a9b7aacf13e25ae47a2ee2f185c14d044ecc3788 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 11:22:56 -0700 Subject: [PATCH 14/21] move onto python 3.12 --- PIQQA.py | 8 ++++---- README.md | 2 +- piqqa-conda-install.txt | 10 +++++----- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 35d604f..4504075 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -865,7 +865,7 @@ def doBoxPlots( if tmpDF.empty: tmpDF = thisDF.copy() else: - tmpDF = tmpDF.append(thisDF, ignore_index=True) + tmpDF = pd.concat([tmpDF, thisDF], ignore_index=True) columnsToKeep = [ "snclq", @@ -976,7 +976,7 @@ def doBoxPlots( print(f" INFO: Generating Boxplots for {channelGroup}") tmpDF = scaledDF[scaledDF["target"].str.endswith(channelGroup)] - grouped = tmpDF.groupby(["station"]) + grouped = tmpDF.groupby("station") filenames = list() for metric in metricList: @@ -1203,7 +1203,7 @@ def doPDFs( # Use the dataframe with the SCALED sample_rms to determine the top/bottom station/target tmpDF = scaledDF[scaledDF["target"].str.endswith(channelGroup)] - grouped = tmpDF.groupby(["snclq"]) + grouped = tmpDF.groupby("snclq") try: df2 = pd.DataFrame( @@ -1451,7 +1451,7 @@ def doSpectrograms( ] # Select the lowest and greatest rms from those stations - grouped = tmpDF.groupby(["snclq"]) + grouped = tmpDF.groupby("snclq") try: df2 = pd.DataFrame( diff --git a/README.md b/README.md index f918e3d..32745c7 100644 --- a/README.md +++ b/README.md @@ -71,7 +71,7 @@ perform installs, updates, or run PIQQA. ``` cd piqqa conda update conda -conda create --name piqqa python=3.9 -c conda-forge --file piqqa-conda-install.txt +conda create --name piqqa python=3.12 -c conda-forge --file piqqa-conda-install.txt conda activate piqqa ``` diff --git a/piqqa-conda-install.txt b/piqqa-conda-install.txt index d7ec76f..29e7fa9 100644 --- a/piqqa-conda-install.txt +++ b/piqqa-conda-install.txt @@ -1,5 +1,5 @@ -pandas=1.3.0 -numpy=1.21.0 -requests= 2.25.1 -matplotlib= 3.4.2 -plotly +pandas=2.1.4 +numpy=1.26.3 +requests=2.31.0 +matplotlib=3.8.2 +plotly=6.0.1 From 31c60c1ccd88bea21816646fdb7e75f3a9175397 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 11:37:11 -0700 Subject: [PATCH 15/21] update IRIS->EarthScope --- PIQQA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PIQQA.py b/PIQQA.py index 4504075..6b9ba5a 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -175,7 +175,7 @@ def doAvailability( actualChannels = sorted(list(set(allMetadataDF["channel"]))) except: quit( - "\nQUITTING: No information found for network and times, is something wrong? Is there data at IRIS for this network? Does this network have channels with the Z component?" + "\nQUITTING: No information found for network and times, is something wrong? Is there data at EarthScope for this network? Does this network have channels with the Z component?" ) ## Loop over the channel groups From 9821ad64985273a8025c679aaed0a5dcfe4a5ef8 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 11:41:17 -0700 Subject: [PATCH 16/21] add check for earthscope OR iris in fedcat return --- reportUtils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/reportUtils.py b/reportUtils.py index 9a3b49a..cc925fc 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -128,9 +128,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): parse_dates=["Earliest", "Latest"], ) tmpDF["staloc"] = f"{s}.{luse}" - availabilityDF = pd.concat( - [availabilityDF, tmpDF], ignore_index=True - ) + availabilityDF = pd.concat([availabilityDF, tmpDF], ignore_index=True) except Exception as e: pass @@ -419,7 +417,7 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve f"starttime={startDate}&endtime={endDate}" f"&format=request&includeoverlaps=false" ) - + print(f" {fedURL}") try: with urllib.request.urlopen(fedURL) as response: html_content = response.read().decode("utf-8") @@ -428,9 +426,11 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve for ln in html_content.split("\n"): if ln.startswith("STATIONSERVICE="): serviceURL = ln.split("=")[1] - if "earthscope" in serviceURL: + if "earthscope" in serviceURL or "iris" in serviceURL: services.append(serviceURL) + print(f"TEMP: services = {services}") + except Exception as e: print( " ERROR: unable to retrieve fed catalog information about where the data lives - %s\n%s " From 6682b2626094a76c8e5b2b4f6d523574592748a7 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 11:42:02 -0700 Subject: [PATCH 17/21] remove temporary logging --- reportUtils.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/reportUtils.py b/reportUtils.py index cc925fc..18d5e28 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -429,8 +429,6 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve if "earthscope" in serviceURL or "iris" in serviceURL: services.append(serviceURL) - print(f"TEMP: services = {services}") - except Exception as e: print( " ERROR: unable to retrieve fed catalog information about where the data lives - %s\n%s " From b1f3b9ab858284d2444a12446a877207c0e75556 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 12:02:53 -0700 Subject: [PATCH 18/21] check availability service for both ph5 and fdsn services --- PIQQA.py | 74 +++++++++++++++++------------------- reportUtils.py | 100 +++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 131 insertions(+), 43 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 6b9ba5a..5368680 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -39,7 +39,6 @@ import numpy as np import shutil import math -import urllib.request def checkAvailability( @@ -128,27 +127,12 @@ def doAvailability( imageDir, tsMetricsExist, ): - # test out whether the availability service is up and running; if not, fall back to + # Test whether the availability service(s) relevant to this network are up (fdsnws-availability + # and ph5ws-availability can go down independently of each other); if not, fall back to # MUSTANG percent_availability extents (no gap-level detail) rather than skipping the plot - availability_exists = True - availablity_test_url = "https://service.earthscope.org/fdsnws/availability/1/" - try: - urllib.request.urlopen(availablity_test_url) - except urllib.error.HTTPError as e: - if e.code == 410: - print( - "WARNING: Availability service is down, falling back to MUSTANG percent_availability for data extents" - ) - else: - print( - f"WARNING: Availability service returned HTTP {e.code}, falling back to MUSTANG percent_availability for data extents" - ) - availability_exists = False - except Exception as e: - print( - f"WARNING: Availability service is down ({e}), falling back to MUSTANG percent_availability for data extents" - ) - availability_exists = False + availability_exists_by_service = reportUtils.checkAvailabilityServicesExist( + network, stations, locations, channels, startDate, endDate + ) print("INFO: Producing availability plot") avFilesDictionary = {} # collects plot filenames @@ -284,6 +268,15 @@ def doAvailability( ) # Only used to get the number of stations, the list isn't actually used otherwise nsta = len(stationList) + # Determine which availability service(s) this channel group's targets actually need + # (fdsnws and/or ph5ws), and only draw the extents/gaps plot if all of them are up. + relevantServices = { + reportUtils.serviceForQuality(s.split(".")[-1]) for s in pctAvDF["snclq"] + } + channelGroupAvailabilityExists = all( + availability_exists_by_service.get(svc, False) for svc in relevantServices + ) + if nsta > nBoxPlotSta: splitPlots = 1 @@ -299,7 +292,7 @@ def doAvailability( topStationsDict[channelGroup] = topStationList - if not availability_exists: + if not channelGroupAvailabilityExists: # Can't draw the extents/gaps plot without the availability service - skip the plot, # but topStationsDict above still lets PDFs/Spectrograms pick stations for this channel group. continue @@ -549,7 +542,7 @@ def doAvailability( ax2.autoscale() ax2.invert_yaxis() ax2.set_yticks(np.arange(nBottom)) - ax2.set_yticklabels(bottfomLabels) + ax2.set_yticklabels(bottomLabels) ax2.set_xlim( [mpl.dates.datestr2num(startDate), mpl.dates.datestr2num(endDate)] ) @@ -588,7 +581,7 @@ def doAvailability( allStationList = pctAvDF["station"].tolist() topStationsDict[channelGroup] = allStationList - if not availability_exists: + if not channelGroupAvailabilityExists: # Can't draw the extents/gaps plot without the availability service - skip the plot, # but topStationsDict above still lets PDFs/Spectrograms pick stations for this channel group. continue @@ -752,7 +745,7 @@ def doAvailability( avFilesDictionary, services_used, availabilityType, - availability_exists, + availability_exists_by_service, ) @@ -1185,7 +1178,7 @@ def doPDFs( network, stations, locations, - availabilityExists, + availabilityExistsByService, ): pdfDictionary = {} @@ -1249,7 +1242,7 @@ def doPDFs( startDate, endDate, tolerance, - availabilityExists, + availabilityExistsByService, ) foundit = False @@ -1322,7 +1315,7 @@ def doPDFs( startDate, endDate, tolerance, - availabilityExists, + availabilityExistsByService, ) foundit = False @@ -1422,7 +1415,7 @@ def doSpectrograms( powerRanges, spectColorPalette, imageDir, - availabilityExists, + availabilityExistsByService, ): spectDictionary = {} # used to track the plots to be used in the final report @@ -1498,7 +1491,7 @@ def doSpectrograms( startDate, endDate, tolerance, - availabilityExists, + availabilityExistsByService, ) foundit = False @@ -1570,7 +1563,7 @@ def doSpectrograms( startDate, endDate, tolerance, - availabilityExists, + availabilityExistsByService, ) foundit = False @@ -1747,7 +1740,6 @@ def doReport( spectColorPalette, powerRanges, userDefinedPowerRange, - availabilityExists, ): print(f"INFO: Writing to file {outfile}") @@ -1797,8 +1789,11 @@ def doReport( # stationServiceLink = f'{service}-station service,' # availabilityServiceLink = f'{service}-availability service,' # stationServiceLink = f"{stationServiceLink[:-1]}" - if availabilityExists: - availabilityServiceLink = f'fdsnws-availability service' + if services: + availabilityServiceLink = " and ".join( + f'{service}-availability service' + for service in services + ) else: availabilityServiceLink = f'mustang service' stationServiceLink = f'fdsnws-station service' @@ -1848,10 +1843,10 @@ def doReport( f"{availabilityIntroText}
      Displaying all stations in the network.
      " ) - if not availabilityExists: + if not services: availabilityIntroText = ( f"{availabilityIntroText}

      " - f"Note: the fdsnws-availability service was unavailable when this report was generated, " + f"Note: the availability service was unavailable when this report was generated, " f"so the data-extent and gap plot(s) below could not be produced. Station selection above was " f"still based on percent-availability metrics from MUSTANG.
      " ) @@ -2828,7 +2823,7 @@ def main(): avFilesDictionary, services, availabilityType, - availabilityExists, + availabilityExistsByService, ] = doAvailability( splitPlots, startDate, @@ -2881,7 +2876,7 @@ def main(): network, stations, locations, - availabilityExists, + availabilityExistsByService, ) # Create Spectrogram Plots @@ -2897,7 +2892,7 @@ def main(): powerRanges, spectColorPalette, imageDir, - availabilityExists, + availabilityExistsByService, ) # Generate Station Map @@ -2945,7 +2940,6 @@ def main(): spectColorPalette, powerRanges, userDefinedPowerRange, - availabilityExists, ) # Zip the report diff --git a/reportUtils.py b/reportUtils.py index 18d5e28..ff19ac0 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -226,15 +226,21 @@ def getDataStartEnd( startDate, endDate, tolerance, - availabilityExists, + availabilityExistsByService, ): """Get the start/end times of the data for a target, via the availability service if it's up, otherwise falling back to MUSTANG extents. Falls back to the report's own startDate/endDate if no start/end can be determined - either way (e.g. percent_availability is also down or has no data).""" + either way (e.g. percent_availability is also down or has no data). + availabilityExistsByService is a dict like {"fdsnws": True, "ph5ws": False} - + the target's own quality code (its last dot-component) picks which entry applies, + since fdsnws-availability and ph5ws-availability can go down independently.""" reportStart = f"{startDate}T00:00:00" reportEnd = f"{endDate}T00:00:00" + serviceKey = serviceForQuality(target.split(".")[-1]) + availabilityExists = availabilityExistsByService.get(serviceKey, False) + if availabilityExists: thisStationAvDF, _ = getAvailability( [target], startDate, endDate, tolerance, "" @@ -370,6 +376,95 @@ def addMetricToDF( return DF +def serviceForQuality(quality): + # PH5/PASSCAL data (quality "D") is served by ph5ws (service.iris.edu), while + # regular FDSN data (quality "M") is served by fdsnws (service.earthscope.org) - + # these can go down independently of each other. + return "ph5ws" if quality == "D" else "fdsnws" + + +def getAvailabilityServiceURLs( + network, stations, locations, channels, startDate, endDate +): + """Determine which availability service(s) - fdsnws and/or ph5ws - actually apply + to this network via FedCatalog, so callers can check/use the correct one(s) instead + of assuming fdsnws. Returns a dict like {"fdsnws": "https://.../availability/1/"}. + """ + chanList = list() + for chan in channels.split(","): + if (len(chan) == 3) and (chan[2] in ["*", "?", ".", "_"]): + chan = f"{chan[0:2]}Z" + elif len(chan) == 2: + chan = f"{chan}Z" + elif chan == "*": + chan = "??Z" + chanList.append(chan) + + fedURL = ( + f"http://service.earthscope.org/irisws/fedcatalog/1/query?" + f"net={network}&sta={stations}&loc={locations}&cha={','.join(chanList)}&" + f"starttime={startDate}&endtime={endDate}" + f"&format=request&includeoverlaps=false" + ) + + serviceURLs = {} + try: + with urllib.request.urlopen(fedURL) as response: + html_content = response.read().decode("utf-8") + for ln in html_content.split("\n"): + if ln.startswith("AVAILABILITYSERVICE="): + url = ln.split("=")[1] + if "ph5ws" in url: + serviceURLs["ph5ws"] = url + elif "fdsnws" in url: + serviceURLs["fdsnws"] = url + except Exception as e: + print( + f" ERROR: unable to retrieve fed catalog information about where the data lives - {fedURL}\n{e}" + ) + # Fall back to checking both - each is probed independently, so this is conservative, not incorrect. + serviceURLs = { + "fdsnws": "https://service.earthscope.org/fdsnws/availability/1/", + "ph5ws": "http://service.iris.edu/ph5ws/availability/1/", + } + + return serviceURLs + + +def checkAvailabilityServicesExist( + network, stations, locations, channels, startDate, endDate +): + """Probe each availability service relevant to this network independently, since + fdsnws-availability and ph5ws-availability can go down independently of each other. + Returns a dict like {"fdsnws": True, "ph5ws": False} covering only the service(s) + actually relevant to this network.""" + serviceURLs = getAvailabilityServiceURLs( + network, stations, locations, channels, startDate, endDate + ) + availabilityExistsByService = {} + for serviceKey, url in serviceURLs.items(): + try: + urllib.request.urlopen(url) + availabilityExistsByService[serviceKey] = True + except urllib.error.HTTPError as e: + if e.code == 410: + print( + f"WARNING: {serviceKey}-availability service is down, falling back to MUSTANG percent_availability for data extents" + ) + else: + print( + f"WARNING: {serviceKey}-availability service returned HTTP {e.code}, falling back to MUSTANG percent_availability for data extents" + ) + availabilityExistsByService[serviceKey] = False + except Exception as e: + print( + f"WARNING: {serviceKey}-availability service is down ({e}), falling back to MUSTANG percent_availability for data extents" + ) + availabilityExistsByService[serviceKey] = False + + return availabilityExistsByService + + def getMetadata(network, stations, locations, channels, startDate, endDate, level): # This one and getStations are almost redundant, except that they return at different # levels. Merge into one function that also takes level as an input? @@ -417,7 +512,6 @@ def getMetadata(network, stations, locations, channels, startDate, endDate, leve f"starttime={startDate}&endtime={endDate}" f"&format=request&includeoverlaps=false" ) - print(f" {fedURL}") try: with urllib.request.urlopen(fedURL) as response: html_content = response.read().decode("utf-8") From 1e9f59adb45d53feff3b5b77912c083f14a3fa47 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 12:36:34 -0700 Subject: [PATCH 19/21] fix typo lowest->highest --- PIQQA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PIQQA.py b/PIQQA.py index 5368680..8e86fc0 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -1540,7 +1540,7 @@ def doSpectrograms( spectDictionary[f"{channelGroup}_smallest"] = sorted(spectFiles) # 2. largest corrected sample_rms station - print(" INFO: Searching for lowest corrected sample_rms station") + print(" INFO: Searching for highest corrected sample_rms station") for ii in reversed(range(len(meds.index))): largestTarget = meds.index[ii] largestNSL = largestTarget[:-3] From dd8caf7ff5bd8fed475717b496d96b6d39202a1a Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 13:48:13 -0700 Subject: [PATCH 20/21] remove future-deprecated params from availability call, create in-house merging for gaps --- reportUtils.py | 52 ++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 6 deletions(-) diff --git a/reportUtils.py b/reportUtils.py index ff19ac0..5b5ef3b 100644 --- a/reportUtils.py +++ b/reportUtils.py @@ -89,6 +89,44 @@ def get_extents_from_mustang(snclqs, startDate, endDate): return availabilityDF +def mergeGapsWithinTolerance(availabilityDF, tolerance): + """Merge adjacent time spans (per network/station/location/channel/quality) that are + separated by a gap no larger than `tolerance` seconds. The availability service used + to do this server-side via the mergegaps parameter; the migrated service no longer + honors it, so this reproduces the same gap-tolerance behavior client-side instead.""" + if availabilityDF.empty: + return availabilityDF + + tolerance = float(tolerance) + groupCols = [ + c + for c in ["network", "station", "location", "channel", "quality"] + if c in availabilityDF.columns + ] + if not groupCols: + return availabilityDF + + mergedRows = [] + for _, group in availabilityDF.groupby(groupCols, dropna=False): + group = group.sort_values("earliest") + current = None + for _, row in group.iterrows(): + if current is None: + current = row.copy() + continue + gapSeconds = (row["earliest"] - current["latest"]).total_seconds() + if gapSeconds <= tolerance: + if row["latest"] > current["latest"]: + current["latest"] = row["latest"] + else: + mergedRows.append(current) + current = row.copy() + if current is not None: + mergedRows.append(current) + + return pd.DataFrame(mergedRows).reset_index(drop=True) + + def getAvailability(snclqs, startDate, endDate, tolerance, avtype): availabilityDF = pd.DataFrame() services = [] @@ -117,8 +155,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): URL = ( f"http://service.earthscope.org/{service}/availability/1/query?format=text&" f"net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&" - f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" - f"mergegaps={tolerance}&includerestricted=true&nodata=404" + f"starttime={startDate}&endtime={endDate}&includerestricted=true&nodata=404" ) try: tmpDF = pd.read_csv( @@ -154,8 +191,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): URL = ( f"http://service.earthscope.org/{service}/availability/1/extent?format=text&" f"net={nets}&sta={stas}&loc={locs}&cha={chans}&quality={qs}&" - f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" - f"includerestricted=true&nodata=404" + f"starttime={startDate}&endtime={endDate}&includerestricted=true&nodata=404" ) try: availabilityDF = pd.read_csv( @@ -192,8 +228,7 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): URL = ( f"http://service.earthscope.org/{service}/availability/1/extent?format=text&" f"net={n}&sta={s}&loc={luse}&cha={c}&quality={q}&" - f"starttime={startDate}&endtime={endDate}&orderby=nslc_time_quality_samplerate&" - f"&includerestricted=true&nodata=404" + f"starttime={startDate}&endtime={endDate}&includerestricted=true&nodata=404" ) try: @@ -214,6 +249,11 @@ def getAvailability(snclqs, startDate, endDate, tolerance, avtype): availabilityDF.rename(columns=lambda x: x.strip().lower(), inplace=True) availabilityDF.rename(columns={"#network": "network"}, inplace=True) + if avtype == "": + # The availability service no longer merges gaps server-side (the mergegaps + # parameter is a no-op on the migrated service), so do it here instead. + availabilityDF = mergeGapsWithinTolerance(availabilityDF, tolerance) + return availabilityDF, services From fed2d93cede932347efba7b964f71455224508f1 Mon Sep 17 00:00:00 2001 From: Laura Keyson Date: Tue, 28 Jul 2026 13:48:29 -0700 Subject: [PATCH 21/21] bypass ts_percent_availability_total if tsmetrics are down --- PIQQA.py | 60 ++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/PIQQA.py b/PIQQA.py index 8e86fc0..4a60760 100755 --- a/PIQQA.py +++ b/PIQQA.py @@ -42,18 +42,20 @@ def checkAvailability( - thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd + thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd, tsMetricsExist ): - thisAvDF = reportUtils.addMetricToDF( - "ts_percent_availability_total", - pd.DataFrame(), - thisNetwork, - thisStation, - thisLocation, - thisChannel, - thisStart, - thisEnd, - ) + thisAvDF = pd.DataFrame() + if tsMetricsExist: + thisAvDF = reportUtils.addMetricToDF( + "ts_percent_availability_total", + pd.DataFrame(), + thisNetwork, + thisStation, + thisLocation, + thisChannel, + thisStart, + thisEnd, + ) if thisAvDF.empty: # need to round out the start and end times for non-ts_percent_availability_total, in case it's very short @@ -1179,6 +1181,7 @@ def doPDFs( stations, locations, availabilityExistsByService, + tsMetricsExist, ): pdfDictionary = {} @@ -1251,7 +1254,13 @@ def doPDFs( ) # get the ts_percent_availatility_total (or percent_availability if ph5) for the time between the data start and end. thisPctAvail = checkAvailability( - thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd + thisNetwork, + thisStation, + thisLocation, + thisChannel, + thisStart, + thisEnd, + tsMetricsExist, ) if thisPctAvail > 75: @@ -1324,7 +1333,13 @@ def doPDFs( ) # get the percent availability between the start and end times of the data thisPctAvail = checkAvailability( - thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd + thisNetwork, + thisStation, + thisLocation, + thisChannel, + thisStart, + thisEnd, + tsMetricsExist, ) if thisPctAvail > 75: @@ -1416,6 +1431,7 @@ def doSpectrograms( spectColorPalette, imageDir, availabilityExistsByService, + tsMetricsExist, ): spectDictionary = {} # used to track the plots to be used in the final report @@ -1500,7 +1516,13 @@ def doSpectrograms( ) # Get the percent availability between the start and end for this target thisPctAvail = checkAvailability( - thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd + thisNetwork, + thisStation, + thisLocation, + thisChannel, + thisStart, + thisEnd, + tsMetricsExist, ) if thisPctAvail > 75: @@ -1572,7 +1594,13 @@ def doSpectrograms( ) # Get the percent availability between the start and end times for this target thisPctAvail = checkAvailability( - thisNetwork, thisStation, thisLocation, thisChannel, thisStart, thisEnd + thisNetwork, + thisStation, + thisLocation, + thisChannel, + thisStart, + thisEnd, + tsMetricsExist, ) if thisPctAvail > 75: @@ -2877,6 +2905,7 @@ def main(): stations, locations, availabilityExistsByService, + tsMetricsExist, ) # Create Spectrogram Plots @@ -2893,6 +2922,7 @@ def main(): spectColorPalette, imageDir, availabilityExistsByService, + tsMetricsExist, ) # Generate Station Map