In [1]:
import time
start_runtime = time.time()

IOOS System Test: Extreme Events Theme: Coastal Inundation

Can we obtain observed current data at stations located within a bounding box?

This notebook is based on IOOS System Test: Inundation

Methodology:

  • Define temporal and spatial bounds of interest, as well as parameters of interest
  • Search for available service endpoints in the NGDC CSW catalog meeting search criteria
  • Search for available OPeNDAP data endpoints
  • Obtain observation data sets from stations within the spatial boundaries (from CO-OPS and NDBC)
  • Extract time series for identified stations
  • Plot time series data, current rose, annual max values per station
  • Plot observation stations on a map

import required libraries


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()


/Users/rpsdev/anaconda/lib/python2.7/site-packages/pandas/io/excel.py:626: UserWarning: Installed openpyxl is not supported at this time. Use >=1.6.1 and <2.0.0.
  .format(openpyxl_compat.start_ver, openpyxl_compat.stop_ver))
Out[2]:
Temporal Bounds - Anything longer than one year kills the CO-OPS service

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))


2004-08-24 14:00 to 2014-08-22 14:00 

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)


26 csw records found
Near-Real Time Surface Ocean Velocity
Near-Real Time Surface Ocean Velocity
Near-Real Time Surface Ocean Velocity
Near-Real Time Surface Ocean Velocity
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 1 km Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
1 km Resolution
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 2 km Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
2 km Resolution
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 6 km Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
6 km Resolution
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 1 km Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
1 km Resolution
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 2 km Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
2 km Resolution
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 6 km Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
6 km Resolution
Surface CUrrents from a Diagnostic model (SCUD): Pacific
Barotropic Tide Model for the Pacific Basin
NOAA.NOS.CO-OPS SOS
National Data Buoy Center SOS
Solar Geophysical Data (SGD) Reports (1955-2009)
HYbrid Coordinate Ocean Model (HYCOM): Global
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 500 m Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
500 m Resolution
Near-Real Time Surface Ocean Velocity, U.S. West Coast, 500 m Resolution (GNOME)
Near-Real Time Surface Ocean Velocity, U.S. West Coast,
500 m Resolution

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[:]))


Total DAP: 11
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/1km/hourly/GNOME
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/1km/hourly/RTV
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/2km/hourly/GNOME
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/2km/hourly/RTV
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/500m/hourly/GNOME
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/500m/hourly/RTV
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/6km/hourly/GNOME
http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USWC/6km/hourly/RTV
http://oos.soest.hawaii.edu/thredds/dodsC/hioos/tide_pac
http://oos.soest.hawaii.edu/thredds/dodsC/pacioos/hycom/global
http://oos.soest.hawaii.edu/thredds/dodsC/pacioos/scud/pac

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))


Total SOS: 2
http://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/SOS?service=SOS&request=GetCapabilities
http://sdf.ndbc.noaa.gov/sos/server.php?service=SOS&request=GetCapabilities&acceptVersions=1.0.0

Update SOS time-date


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')
Get list of stations - we get a list of the available stations from NOAA and COOPS

Initialize Station Data List


In [10]:
st_list = {}

Get CO-OPS Station Data


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))


2004-08-24 14:00:00+00:00:2014-08-22 14:00:00+00:00
1023

gets a list of the active stations from coops


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)


http://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/SOS?service=SOS&request=GetObservation&version=1.0.0&observedProperty=Currents&bin=1&offering=urn:ioos:network:NOAA.NOS.CO-OPS:CurrentsActive&featureOfInterest=BBOX:-130,38,-121,50&responseFormat=text/csv
Date: 2004-08-24T14:00:00Z to 2014-08-22T14:00:00Z
Lat/Lon Box: -130,38,-121,50

COOPS Station Information


In [13]:
st_list = processStationInfo(obs_loc_df, st_list, "coops")


urn:ioos:station:NOAA.NOS.CO-OPS:cp0101
urn:ioos:station:NOAA.NOS.CO-OPS:hb0401
urn:ioos:station:NOAA.NOS.CO-OPS:s06010
number of stations in bbox 3

In [14]:
st_list


Out[14]:
{'urn:ioos:station:NOAA.NOS.CO-OPS:cp0101': {'lat': 48.8628,
  'lon': -122.761,
  'source': 'coops'},
 'urn:ioos:station:NOAA.NOS.CO-OPS:hb0401': {'lat': 40.777500000000003,
  'lon': -124.1966,
  'source': 'coops'},
 'urn:ioos:station:NOAA.NOS.CO-OPS:s06010': {'lat': 38.034599999999998,
  'lon': -122.12520000000001,
  'source': 'coops'}}

Get NDBC Station Data


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))


2004-08-24 14:00:00+00:00:2014-08-22 14:00:00+00:00
891

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)


Date: 2004-08-24T14:00:00Z to 2014-08-22T14:00:00Z
Lat/Lon Box: -130,38,-121,50
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:-130,38,-121,50&observedproperty=Currents&responseformat=text/csv&

NDBC Station information


In [17]:
st_list = processStationInfo(obs_loc_df, st_list, "ndbc")
st_list


urn:ioos:station:wmo:46013
urn:ioos:station:wmo:46015
urn:ioos:station:wmo:46022
urn:ioos:station:wmo:46029
urn:ioos:station:wmo:46087
urn:ioos:station:wmo:46088
urn:ioos:station:wmo:46089
number of stations in bbox 10
Out[17]:
{'urn:ioos:station:NOAA.NOS.CO-OPS:cp0101': {'lat': 48.8628,
  'lon': -122.761,
  'source': 'coops'},
 'urn:ioos:station:NOAA.NOS.CO-OPS:hb0401': {'lat': 40.777500000000003,
  'lon': -124.1966,
  'source': 'coops'},
 'urn:ioos:station:NOAA.NOS.CO-OPS:s06010': {'lat': 38.034599999999998,
  'lon': -122.12520000000001,
  'source': 'coops'},
 'urn:ioos:station:wmo:46013': {'lat': 38.229999999999997,
  'lon': -123.31999999999999,
  'source': 'ndbc'},
 'urn:ioos:station:wmo:46015': {'lat': 42.759999999999998,
  'lon': -124.83,
  'source': 'ndbc'},
 'urn:ioos:station:wmo:46022': {'lat': 40.780000000000001,
  'lon': -124.54000000000001,
  'source': 'ndbc'},
 'urn:ioos:station:wmo:46029': {'lat': 46.140000000000001,
  'lon': -124.51000000000001,
  'source': 'ndbc'},
 'urn:ioos:station:wmo:46087': {'lat': 49.890000000000001,
  'lon': -124.73,
  'source': 'ndbc'},
 'urn:ioos:station:wmo:46088': {'lat': 48.329999999999998,
  'lon': -123.17,
  'source': 'ndbc'},
 'urn:ioos:station:wmo:46089': {'lat': 45.909999999999997,
  'lon': -125.76000000000001,
  'source': 'ndbc'}}

In [18]:
print(st_list[st_list.keys()[0]]['lat'])
print(st_list[st_list.keys()[0]]['lon'])


38.0346
-122.1252

The function only support who date time differences

Large Temporal Requests Need To Be Broken Down - When requesting a large temporal range outside the SOS limit, the sos request needs to be broken down. See issues in [ioos](https://github.com/ioos/system-test/issues/81), [ioos](https://github.com/ioos/system-test/issues/101), [ioos](https://github.com/ioos/system-test/issues/116) and [pyoos](https://github.com/ioos/pyoos/issues/35). Unfortunately currents is not available via DAP ([ioos](https://github.com/ioos/system-test/issues/116))
Large Temporal Requests Need To Be Broken Down - Obtaining long time series from COOPS via SOS is not ideal and the opendap links are not available, so we use the tides and currents api to get the currents in json format. The api response provides in default bin, unless a bin is specified (i.e bin=1)
Pyoos - Should be able to use the collector, but does not work?
Use NDBC DAP endpoints to get time-series data - The DAP server for currents is available for NDBC data, we use that to get long time series data.
Progress Information For Large Requests - Shows the user a progress bar for each stations as its processed. Click [here]('http://www.tidesandcurrents.noaa.gov/cdata/StationList?type=Current+Data&filter=active') to show more information on the CO-OPS locations
Processing long time series - The CO-OPS Server responds really slow (> 30 secs, for what should be a 5 sec request) to multiple requests, so getting long time series data is almost impossible.

get CO-OPS station data


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


start & end dates: 2004-08-24 14:00:00, 2014-08-22 14:00:00

[coops]: urn:ioos:station:NOAA.NOS.CO-OPS:s06010
[coops]: urn:ioos:station:NOAA.NOS.CO-OPS:cp0101
[ndbc]: urn:ioos:station:wmo:46022
[coops]: urn:ioos:station:NOAA.NOS.CO-OPS:hb0401
[ndbc]: urn:ioos:station:wmo:46087
[ndbc]: urn:ioos:station:wmo:46015
[ndbc]: urn:ioos:station:wmo:46029
[ndbc]: urn:ioos:station:wmo:46013
[ndbc]: urn:ioos:station:wmo:46089
[ndbc]: urn:ioos:station:wmo:46088

In [20]:
# Check theres data in there.
st_list[st_list.keys()[2]]


Out[20]:
{'hasObsData': True,
 'lat': 40.780000000000001,
 'lon': -124.54000000000001,
 'obsData':                      sea_water_speed (cm/s)  \
 2010-09-01 00:00:00                       7   
 2010-09-01 01:00:00                       7   
 2010-09-01 02:00:00                       4   
 2010-09-01 03:00:00                       8   
 2010-09-01 04:00:00                      11   
 2010-09-01 05:00:00                      24   
 2010-09-01 06:00:00                      42   
 2010-09-01 07:00:00                      51   
 2010-09-01 08:00:00                      51   
 2010-09-01 09:00:00                      36   
 2010-09-01 10:00:00                      40   
 2010-09-01 11:00:00                      37   
 2010-09-01 12:00:00                      35   
 2010-09-01 13:00:00                      37   
 2010-09-01 14:00:00                      37   
 2010-09-01 15:00:00                      34   
 2010-09-01 16:00:00                      19   
 2010-09-01 17:00:00                      19   
 2010-09-01 18:00:00                       9   
 2010-09-01 19:00:00                      17   
 2010-09-01 20:00:00                      24   
 2010-09-01 21:00:00                      27   
 2010-09-01 22:00:00                      27   
 2010-09-01 23:00:00                      24   
 2010-09-02 00:00:00                      20   
 2010-09-02 01:00:00                      22   
 2010-09-02 02:00:00                      22   
 2010-09-02 03:00:00                      26   
 2010-09-02 04:00:00                      35   
 2010-09-02 05:00:00                      37   
 ...                                     ...   
 2011-01-02 04:00:00                     NaN   
 2011-01-03 00:00:00                     NaN   
 2011-01-03 03:00:00                     NaN   
 2011-01-03 10:00:00                     NaN   
 2011-01-04 01:00:00                     NaN   
 2011-01-04 10:00:00                     NaN   
 2011-01-04 12:00:00                     NaN   
 2011-01-05 00:00:00                     NaN   
 2011-01-05 02:00:00                     NaN   
 2011-01-06 13:00:00                     NaN   
 2011-01-06 15:00:00                     NaN   
 2011-01-06 19:00:00                     NaN   
 2011-01-07 02:00:00                     NaN   
 2011-01-07 15:00:00                     NaN   
 2011-01-07 17:00:00                     NaN   
 2011-01-08 11:00:00                     NaN   
 2011-01-08 16:00:00                     NaN   
 2011-01-10 18:00:00                     NaN   
 2011-01-11 10:00:00                     NaN   
 2011-01-11 16:00:00                     NaN   
 2011-01-11 20:00:00                     NaN   
 2011-01-12 19:00:00                     NaN   
 2011-01-13 01:00:00                     NaN   
 2011-01-13 07:00:00                     NaN   
 2011-01-13 11:00:00                     NaN   
 2011-01-14 09:00:00                     NaN   
 2011-01-14 22:00:00                     NaN   
 2011-01-16 05:00:00                     NaN   
 2011-01-16 07:00:00                     NaN   
 2011-01-23 13:00:00                     NaN   
 
                      direction_of_sea_water_velocity (degree)  
 2010-09-01 00:00:00                                       210  
 2010-09-01 01:00:00                                       177  
 2010-09-01 02:00:00                                       190  
 2010-09-01 03:00:00                                       158  
 2010-09-01 04:00:00                                       142  
 2010-09-01 05:00:00                                       125  
 2010-09-01 06:00:00                                       122  
 2010-09-01 07:00:00                                       135  
 2010-09-01 08:00:00                                       154  
 2010-09-01 09:00:00                                       174  
 2010-09-01 10:00:00                                       184  
 2010-09-01 11:00:00                                       180  
 2010-09-01 12:00:00                                       175  
 2010-09-01 13:00:00                                       190  
 2010-09-01 14:00:00                                       217  
 2010-09-01 15:00:00                                       222  
 2010-09-01 16:00:00                                       239  
 2010-09-01 17:00:00                                       237  
 2010-09-01 18:00:00                                       199  
 2010-09-01 19:00:00                                       153  
 2010-09-01 20:00:00                                       178  
 2010-09-01 21:00:00                                       183  
 2010-09-01 22:00:00                                       188  
 2010-09-01 23:00:00                                       188  
 2010-09-02 00:00:00                                       195  
 2010-09-02 01:00:00                                       192  
 2010-09-02 02:00:00                                       193  
 2010-09-02 03:00:00                                       208  
 2010-09-02 04:00:00                                       201  
 2010-09-02 05:00:00                                       192  
 ...                                                       ...  
 2011-01-02 04:00:00                                        16  
 2011-01-03 00:00:00                                       224  
 2011-01-03 03:00:00                                        16  
 2011-01-03 10:00:00                                        16  
 2011-01-04 01:00:00                                       227  
 2011-01-04 10:00:00                                        16  
 2011-01-04 12:00:00                                        16  
 2011-01-05 00:00:00                                        16  
 2011-01-05 02:00:00                                        16  
 2011-01-06 13:00:00                                        16  
 2011-01-06 15:00:00                                        16  
 2011-01-06 19:00:00                                        16  
 2011-01-07 02:00:00                                        16  
 2011-01-07 15:00:00                                        16  
 2011-01-07 17:00:00                                        16  
 2011-01-08 11:00:00                                        16  
 2011-01-08 16:00:00                                        16  
 2011-01-10 18:00:00                                        16  
 2011-01-11 10:00:00                                        16  
 2011-01-11 16:00:00                                        16  
 2011-01-11 20:00:00                                        16  
 2011-01-12 19:00:00                                        16  
 2011-01-13 01:00:00                                        16  
 2011-01-13 07:00:00                                        16  
 2011-01-13 11:00:00                                        16  
 2011-01-14 09:00:00                                        16  
 2011-01-14 22:00:00                                        16  
 2011-01-16 05:00:00                                        16  
 2011-01-16 07:00:00                                        16  
 2011-01-23 13:00:00                                        16  
 
 [2376 rows x 2 columns],
 'source': 'ndbc'}

Plot the pandas data frames for the stations

Station Data Plot - There might be an issue with some of the NDBC station data...

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


num rows: 2376
num rows: 62073
num rows: 4129
num rows: 10001
num rows: 2153
num rows: 11071
num rows: 97899

Find the min and max data values

Station Data Plot - Some stations might not plot due to the data.

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')


Produce Interactive Map


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]:
"); map.addLayer(marker_3) var marker_4_icon = L.AwesomeMarkers.icon({ icon: 'remove',markerColor: 'red',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_4 = L.marker([40.7775, -124.1966], {'icon':marker_4_icon} ); marker_4.bindPopup("Obs Location:
urn:ioos:station:NOAA.NOS.CO-OPS:hb0401"); map.addLayer(marker_4) var marker_5_icon = L.AwesomeMarkers.icon({ icon: 'ok',markerColor: 'green',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_5 = L.marker([49.89, -124.73], {'icon':marker_5_icon} ); marker_5.bindPopup("Obs Location:
urn:ioos:station:wmo:46087
"); map.addLayer(marker_5) var marker_6_icon = L.AwesomeMarkers.icon({ icon: 'ok',markerColor: 'green',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_6 = L.marker([42.76, -124.83], {'icon':marker_6_icon} ); marker_6.bindPopup("Obs Location:
urn:ioos:station:wmo:46015
"); map.addLayer(marker_6) var marker_7_icon = L.AwesomeMarkers.icon({ icon: 'ok',markerColor: 'green',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_7 = L.marker([46.14, -124.51], {'icon':marker_7_icon} ); marker_7.bindPopup("Obs Location:
urn:ioos:station:wmo:46029
"); map.addLayer(marker_7) var marker_8_icon = L.AwesomeMarkers.icon({ icon: 'ok',markerColor: 'green',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_8 = L.marker([38.23, -123.32], {'icon':marker_8_icon} ); marker_8.bindPopup("Obs Location:
urn:ioos:station:wmo:46013
"); map.addLayer(marker_8) var marker_9_icon = L.AwesomeMarkers.icon({ icon: 'ok',markerColor: 'green',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_9 = L.marker([45.91, -125.76], {'icon':marker_9_icon} ); marker_9.bindPopup("Obs Location:
urn:ioos:station:wmo:46089
"); map.addLayer(marker_9) var marker_10_icon = L.AwesomeMarkers.icon({ icon: 'ok',markerColor: 'green',prefix: 'glyphicon',extraClasses: 'fa-rotate-0'}); var marker_10 = L.marker([48.33, -123.17], {'icon':marker_10_icon} ); marker_10.bindPopup("Obs Location:
urn:ioos:station:wmo:46088
"); map.addLayer(marker_10) var latLngs = [ [38, -130], [38, -121], [50, -121], [50, -130], [38, -130], ]; var line_1 = L.polyline(latLngs,{ color: '#FF0000', weight: 5, }); map.addLayer(line_1); " style="width: 100%; height: 500px; border: none">

In [25]:
elapsed = time.time() - start_runtime
print('{:.2f} minutes'.format(elapsed / 60.))


1.87 minutes