# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <markdowncell>

# # IOOS System Test: [Extreme Events Theme]
# (https://github.com/ioos/system-test/wiki/Development-of-Test-Themes#
# wiki-theme-2-extreme-events):  Inundation
# 
# ### Compare modeled water levels with observations for a specified bounding box
# and time period using IOOS recommended service standards for catalog search
# (CSW) and data retrieval (OPeNDAP & SOS).
# 
# * Query CSW to find datasets that match criteria
# * Extract OPeNDAP data endpoints from model datasets and SOS endpoints from
#   observational datasets
# * OPeNDAP model datasets will be granules
# * SOS endpoints may be datasets (from ncSOS) or collections of datasets (from
#   NDBC, CO-OPS SOS servers)
# * Filter SOS services to obtain datasets
# * Extract data from SOS datasets
# * Extract data from model datasets at locations of observations
# * Compare time series data on same vertical datum

# <codecell>

# Standard Library.
from warnings import warn
from datetime import datetime, timedelta

# Scientific stack.
import iris
iris.FUTURE.netcdf_promote = True

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.transforms import offset_copy
from cartopy.io.img_tiles import MapQuestOpenAerial
from pandas import DataFrame, date_range, read_csv, concat

from iris.exceptions import CoordinateNotFoundError, ConstraintMismatchError

# Custom IOOS/ASA modules (available at PyPI).
from owslib import fes
from owslib.csw import CatalogueServiceWeb
from pyoos.collectors.coops.coops_sos import CoopsSos

# Local imports
from utilities import name_list, sos_name
from utilities import (fes_date_filter, coops2df, find_timevar, find_ij, nearxy,
                       service_urls, mod_df)

# <markdowncell>

# ### Specify a time range and bounding box of interest:

# <codecell>

dates = {'Hurricane sandy':
         [datetime(2012, 10, 26), datetime(2012, 11, 2)],
         '2014 Feb 10-15 Storm':
         [datetime(2014, 2, 10), datetime(2014, 2, 15)],
         '2014 Recent': [datetime(2014, 3, 8), datetime(2014, 3, 11)],
         '2011': [datetime(2013, 4, 20), datetime(2013, 4, 24)]}

jd_now = datetime.utcnow()
jd_start,  jd_stop = jd_now - timedelta(days=3), jd_now + timedelta(days=3)

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

# <codecell>

# 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]}

box = area['New York harbor region']

# <markdowncell>

# ### Search CSW for datasets of interest

# <codecell>

if False:
    from IPython.core.display import HTML
    url = 'http://www.ngdc.noaa.gov/geoportal/'
    HTML('<iframe src=%s width=950 height=400></iframe>' % url)

# <codecell>

# Connect to CSW, explore it's properties.
CSW = {'NGDC Geoportal':
       'http://www.ngdc.noaa.gov/geoportal/csw',
       'USGS WHSC Geoportal':
       'http://geoport.whoi.edu/geoportal/csw',
       'NODC Geoportal: granule level':
       'http://www.nodc.noaa.gov/geoportal/csw',
       'NODC Geoportal: collection level':
       'http://data.nodc.noaa.gov/geoportal/csw',
       'NRCAN CUSTOM':
       'http://geodiscover.cgdi.ca/wes/serviceManagerCSW/csw',
       'USGS Woods Hole GI_CAT':
       'http://geoport.whoi.edu/gi-cat/services/cswiso',
       'USGS CIDA Geonetwork':
       'http://cida.usgs.gov/gdp/geonetwork/srv/en/csw',
       'USGS Coastal and Marine Program':
       'http://cmgds.marine.usgs.gov/geonetwork/srv/en/csw',
       'USGS Woods Hole Geoportal':
       'http://geoport.whoi.edu/geoportal/csw',
       'CKAN testing site for new Data.gov':
       'http://geo.gov.ckan.org/csw',
       'EPA':
       'https://edg.epa.gov/metadata/csw',
       'CWIC':
       'http://cwic.csiss.gmu.edu/cwicv1/discovery'}

endpoint = CSW['NGDC Geoportal']
csw = CatalogueServiceWeb(endpoint, timeout=60)
csw.version

# <codecell>

csw.get_operation_by_name('GetRecords').constraints

# <markdowncell>

# ### Convert User Input into FES filters.

# <codecell>

start, stop = fes_date_filter(start_date, stop_date)
bbox = fes.BBox(box)

# <codecell>

or_filt = fes.Or([fes.PropertyIsLike(propertyname='apiso:AnyText',
                                     literal=('*%s*' % val),
                                     escapeChar='\\',
                                     wildCard='*',
                                     singleChar='?') for val in name_list])

# <markdowncell>

# ROMS model output often has Averages and History files.  The Averages files are usually averaged over a tidal cycle or more, while the History files are snapshots at that time instant.  We are not interested in averaged data for this test, so in the cell below we remove any Averages files here by removing any datasets that have the term "Averages" in the metadata text.  A better approach would be to look at the `cell_methods` attributes propagated through to some term in the ISO metadata, but this is not implemented yet, as far as I know

# <codecell>

val = 'Averages'
not_filt = fes.Not([fes.PropertyIsLike(propertyname='apiso:AnyText',
                                       literal=('*%s*' % val),
                                       escapeChar='\\',
                                       wildCard='*',
                                       singleChar='?')])

filter_list = [fes.And([bbox, start, stop, or_filt, not_filt])]

# <markdowncell>

# Try request using multiple filters "and" syntax: [[filter1,filter2]].

# <codecell>

csw.getrecords2(constraints=filter_list, maxrecords=1000, esn='full')
print(len(csw.records.keys()))

# <markdowncell>

# Now print out some titles

# <codecell>

for rec, item in csw.records.items():
    print(item.title)

# <markdowncell>

# Print out all the OPeNDAP Data URL endpoints

# <codecell>

dap_urls = service_urls(csw.records,
                        service='odp:url')
print("\n".join(dap_urls))

# <markdowncell>

# Print out all the SOS Data URL endpoints

# <codecell>

sos_urls = service_urls(csw.records,
                        service='sos:url')