In [1]:
import time
start_runtime = time.time()
IOOS System Test: Extreme Events Theme: Coastal Inundation
Methodology:
In [2]:
import os
import os.path
from datetime import datetime, timedelta
import uuid
import folium
import matplotlib.pyplot as plt
from owslib.csw import CatalogueServiceWeb
from owslib import fes
import numpy as np
from pandas import read_csv
from pyoos.collectors.ndbc.ndbc_sos import NdbcSos
from pyoos.collectors.coops.coops_sos import CoopsSos
from utilities import (fes_date_filter, service_urls, get_coordinates,
inline_map, css_styles, processStationInfo,
get_ncfiles_catalog, new_axes, set_legend)
css_styles()
Out[2]:
In [3]:
bounding_box_type = "box"
# Bounding Box [lon_min, lat_min, lon_max, lat_max]
area = {'Hawaii': [-160.0, 18.0, -154., 23.0],
'Gulf of Maine': [-72.0, 41.0, -69.0, 43.0],
'New York harbor region': [-75., 39., -71., 41.5],
'Puerto Rico': [-71, 14, -60, 24],
'East Coast': [-77, 34, -70, 40],
'North West': [-130, 38, -121, 50]}
bounding_box = area['North West']
# Temporal range.
jd_now = datetime.utcnow()
jd_start, jd_stop = jd_now - timedelta(days=(365*10)), jd_now
start_date = jd_start.strftime('%Y-%m-%d %H:00')
stop_date = jd_stop.strftime('%Y-%m-%d %H:00')
jd_start = datetime.strptime(start_date, '%Y-%m-%d %H:%M')
jd_stop = datetime.strptime(stop_date, '%Y-%m-%d %H:%M')
print('%s to %s ' % (start_date, stop_date))
In [4]:
# Put the names in a dict for ease of access.
data_dict = {}
sos_name = 'Currents'
data_dict['currents'] = {"names": ['currents',
'surface_eastward_sea_water_velocity',
'*surface_eastward_sea_water_velocity*'],
"sos_name": ['currents']}
CSW Search
In [5]:
endpoint = 'http://www.ngdc.noaa.gov/geoportal/csw' # NGDC Geoportal.
csw = CatalogueServiceWeb(endpoint, timeout=60)
Search
In [6]:
# Convert User Input into FES filters.
start, stop = fes_date_filter(start_date, stop_date)
bbox = fes.BBox(bounding_box)
# Use the search name to create search filter.
kw = dict(propertyname='apiso:AnyText', escapeChar='\\',
wildCard='*', singleChar='?')
or_filt = fes.Or([fes.PropertyIsLike(literal=('*%s*' % val), **kw) for
val in data_dict['currents']['names']])
val = 'Averages'
not_filt = fes.Not([fes.PropertyIsLike(literal=('*%s*' % val), **kw)])
filter_list = [fes.And([bbox, start, stop, or_filt, not_filt])]
# Connect to CSW, explore it's properties
# try request using multiple filters "and" syntax: [[filter1, filter2]]
csw.getrecords2(constraints=filter_list, maxrecords=1000, esn='full')
print("%s csw records found" % len(csw.records))
for rec, item in csw.records.items():
print(item.title)
DAP
In [7]:
dap_urls = service_urls(csw.records)
# Remove duplicates and organize.
dap_urls = sorted(set(dap_urls))
print("Total DAP: %s" % len(dap_urls))
# Print the first 5:
print("\n".join(dap_urls[:]))
Get SOS links, NDBC is not available so add it...
In [8]:
sos_urls = service_urls(csw.records, service='sos:url')
# Remove duplicates and organize.
sos_urls = sorted(set(sos_urls))
print("Total SOS: %s" % len(sos_urls))
print("\n".join(sos_urls))
In [9]:
start_time = datetime.strptime(start_date, '%Y-%m-%d %H:%M')
end_time = datetime.strptime(stop_date, '%Y-%m-%d %H:%M')
iso_start = start_time.strftime('%Y-%m-%dT%H:%M:%SZ')
iso_end = end_time.strftime('%Y-%m-%dT%H:%M:%SZ')
In [10]:
st_list = {}
In [11]:
coops_collector = CoopsSos()
coops_collector.start_time = start_time
coops_collector.end_time = end_time
coops_collector.variables = data_dict["currents"]["sos_name"]
coops_collector.server.identification.title
ofrs = coops_collector.server.offerings
print("%s:%s" % (coops_collector.start_time, coops_collector.end_time))
print(len(ofrs))
In [12]:
box_str = ','.join(str(e) for e in bounding_box)
url = (('http://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/SOS?'
'service=SOS&request=GetObservation&version=1.0.0&'
'observedProperty=%s&bin=1&'
'offering=urn:ioos:network:NOAA.NOS.CO-OPS:CurrentsActive&'
'featureOfInterest=BBOX:%s&responseFormat=text/csv') %
(sos_name, box_str))
obs_loc_df = read_csv(url)
print(url)
print("Date: %s to %s" % (iso_start, iso_end))
print("Lat/Lon Box: %s" % box_str)
In [13]:
st_list = processStationInfo(obs_loc_df, st_list, "coops")
In [14]:
st_list
Out[14]:
In [15]:
ndbc_collector = NdbcSos()
ndbc_collector.start_time = start_time
ndbc_collector.end_time = end_time
ndbc_collector.variables = data_dict["currents"]["sos_name"]
ndbc_collector.server.identification.title
print("%s:%s" % (ndbc_collector.start_time, ndbc_collector.end_time))
ofrs = ndbc_collector.server.offerings
print(len(ofrs))
In [16]:
print("Date: %s to %s" % (iso_start, iso_end))
box_str = ','.join(str(e) for e in bounding_box)
print("Lat/Lon Box: %s" % box_str)
url = (('http://sdf.ndbc.noaa.gov/sos/server.php?'
'request=GetObservation&service=SOS&'
'version=1.0.0&'
'offering=urn:ioos:network:noaa.nws.ndbc:all&'
'featureofinterest=BBOX:%s&'
'observedproperty=%s&'
'responseformat=text/csv&') % (box_str, sos_name))
print(url)
obs_loc_df = read_csv(url)
In [17]:
st_list = processStationInfo(obs_loc_df, st_list, "ndbc")
st_list
Out[17]:
In [18]:
print(st_list[st_list.keys()[0]]['lat'])
print(st_list[st_list.keys()[0]]['lon'])
In [19]:
# Used to define the number of days allowable by the service.
coops_point_max_days = ndbc_point_max_days = 30
print("start & end dates: %s, %s\n" % (jd_start, jd_stop))
for station_index in st_list.keys():
# Set it so we can use it later.
st = station_index.split(":")[-1]
print('[%s]: %s' % (st_list[station_index]['source'], station_index))
divid = str(uuid.uuid4())
if st_list[station_index]['source'] == 'coops':
# Coops fails for large requests.
master_df = []
elif st_list[station_index]['source'] == 'ndbc':
# Use the dap catalog to get the data.
master_df = get_ncfiles_catalog(station_index, jd_start, jd_stop)
if len(master_df) > 0:
st_list[station_index]['hasObsData'] = True
st_list[station_index]['obsData'] = master_df
In [20]:
# Check theres data in there.
st_list[st_list.keys()[2]]
Out[20]:
In [21]:
for station_index in st_list.keys():
df = st_list[station_index]['obsData']
if len(df) > 1:
st_list[station_index]['hasObsData'] = True
print("num rows: %s" % len(df))
fig = plt.figure(figsize=(18, 3))
plt.scatter(df.index, df['sea_water_speed (cm/s)'])
fig.suptitle('Station:'+station_index, fontsize=20)
plt.xlabel('Date', fontsize=18)
plt.ylabel('sea_water_speed (cm/s)', fontsize=16)
else:
st_list[station_index]['hasObsData'] = False
In [22]:
# Build current roses.
filelist = [f for f in os.listdir("./images") if f.endswith(".png")]
for f in filelist:
os.remove("./images/{}".format(f))
station_min_max = {}
for station_index in st_list.keys():
all_spd_data = {}
all_dir_data = {}
all_time_spd = []
all_time_dir = []
df = st_list[station_index]['obsData']
if len(df) > 1:
try:
spd_data = df['sea_water_speed (cm/s)'].values
spd_data = np.array(spd_data)
dir_data = df['direction_of_sea_water_velocity (degree)'].values
dir_data = np.array(dir_data)
time_data = df.index.tolist()
time_data = np.array(time_data)
# NOTE: This data cleanup can a vectorized function.
for idx in range(0, len(spd_data)):
if spd_data[idx] > 998:
continue
elif np.isnan(spd_data[idx]):
continue
elif dir_data[idx] == 0:
continue
else:
dt_year = time_data[idx].year
dt_year = str(dt_year)
if dt_year not in all_spd_data.keys():
all_spd_data[dt_year] = []
all_dir_data[dt_year] = []
# Convert to knots.
knot_val = (spd_data[idx] * 0.0194384449)
knot_val = "%.4f" % knot_val
knot_val = float(knot_val)
all_spd_data[dt_year].append(knot_val)
all_dir_data[dt_year].append(dir_data[idx])
all_time_spd.append(knot_val)
all_time_dir.append(dir_data[idx])
all_time_spd = np.array(all_time_spd, dtype=np.float)
all_time_dir = np.array(all_time_dir, dtype=np.float)
station_min_max[station_index] = {}
for year in all_spd_data.keys():
year_spd = np.array(all_spd_data[year])
year_dir = np.array(all_dir_data[year])
station_min_max[station_index][year] = {}
station_min_max[station_index][year]['pts'] = len(year_spd)
min_spd, max_spd = np.min(year_spd), np.max(year_spd)
station_min_max[station_index][year]['spd_min'] = min_spd
station_min_max[station_index][year]['spd_max'] = max_spd
dir_min, dir_max = np.argmin(year_spd), np.argmax(year_spd)
yr_dir_min, yr_dir_max = year_dir[dir_min], year_dir[dir_max]
station_min_max[station_index][year]['dir_at_min'] = yr_dir_min
station_min_max[station_index][year]['dir_at_max'] = yr_dir_max
try:
# A stacked histogram with normed
# (displayed in percent) results.
ax = new_axes()
ax.set_title(station_index.split(":")[-1] +
" stacked histogram with normed (displayed in %)"
"\nresults (spd in knots), All Time.")
ax.bar(all_time_dir, all_time_spd, normed=True,
opening=0.8, edgecolor='white')
set_legend(ax)
fig = plt.gcf()
fig.set_size_inches(8, 8)
fname = './images/%s.png' % station_index.split(":")[-1]
fig.savefig(fname, dpi=100)
except Exception as e:
print("Error when plotting %s" % e)
pass
except Exception as e: # Be specific here!
print("Error: %s" % e)
pass
In [34]:
# Plot the min and max from each station.
fields = ['spd_']
for idx in range(0, len(fields)):
d_field = fields[idx]
fig, ax = plt.subplots(1, 1, figsize=(18, 5))
for st in station_min_max:
x, y_min, y_max = [], [], []
for year in station_min_max[st]:
x.append(year)
y_max.append(station_min_max[st][year][d_field+'max'])
marker_size = station_min_max[st][year]['pts'] / 80
marker_size += 20
station_label = st.split(":")[-1]
ax.scatter(np.array(x), np.array(y_max),
label=station_label, s=marker_size,
c=np.random.rand(3, 1), marker="o")
ax.set_xlim([2000, 2015])
ax.set_title("Yearly Max Speed Per Station, Marker Scaled Per "
"Annual Pts (bigger = more pts per year)")
ax.set_ylabel("speed (knots)")
ax.set_xlabel("Year")
ax.legend(loc='upper left')
In [24]:
station = st_list[st_list.keys()[0]]
m = folium.Map(location=[station["lat"], station["lon"]], zoom_start=4)
m.line(get_coordinates(bounding_box, bounding_box_type),
line_color='#FF0000', line_weight=5)
# Plot the obs station.
for st in st_list:
hasObs = st_list[st]['hasObsData']
if hasObs:
fname = './images/%s.png' % st.split(":")[-1]
if os.path.isfile(fname):
popup = ('Obs Location:<br>%s<br><img border=120 src="'
'./images/%s.png" width="242" height="242">' %
(st, st.split(":")[-1]))
m.simple_marker([st_list[st]["lat"], st_list[st]["lon"]],
popup=popup,
marker_color="green",
marker_icon="ok")
else:
popup = 'Obs Location:<br>%s' % st
m.simple_marker([st_list[st]["lat"], st_list[st]["lon"]],
popup=popup,
marker_color="green",
marker_icon="ok")
else:
popup = 'Obs Location:<br>%s' % st
m.simple_marker([st_list[st]["lat"], st_list[st]["lon"]],
popup=popup,
marker_color="red",
marker_icon="remove")
inline_map(m)
Out[24]:
In [25]:
elapsed = time.time() - start_runtime
print('{:.2f} minutes'.format(elapsed / 60.))