In [1]:

    

from pylab import *
import datetime as dt
import iris


    

In [2]:

    

#box=[-74.4751, 40.3890, -73.7432, 40.9397]
box=[-76.4751, 38.3890, -71.7432, 42.9397]

#box=[-180, -90, 180, 90]

# specific specific times (UTC) ...

# hurricane sandy
jd_start = dt.datetime(2012,10,26)
jd_stop = dt.datetime(2012,11,2)

# 2014 feb 10-15 storm
jd_start = dt.datetime(2014,2,10)
jd_stop = dt.datetime(2014,2,15)

# 2014 recent
jd_start = dt.datetime(2014,3,8)
jd_stop = dt.datetime(2014,3,11)

# 2011 
#jd_start = dt.datetime(2013,4,20)
#jd_stop = dt.datetime(2013,4,24)

# ... or relative to now
jd_start = dt.datetime.utcnow()- dt.timedelta(days=3)
jd_stop = dt.datetime.utcnow() + dt.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 = dt.datetime.strptime(start_date,'%Y-%m-%d %H:%M')
jd_stop = dt.datetime.strptime(stop_date,'%Y-%m-%d %H:%M')

print start_date,'to',stop_date

sos_name = 'water_surface_height_above_reference_datum'


    

In [3]:

    

url='http://sos.maracoos.org/stable/dodsC/stationHRMARPH-agg.ncml'
url='http://geoport-dev.whoi.edu/thredds/dodsC/usgs/data2/rsignell/foo.ncml'
#url='http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/DBOFS/fmrc/Aggregated_7_day_DBOFS_Fields_Forecast_best.ncd'


    

In [4]:

    

std_name_list=['water_surface_height_above_reference_datum',
    'sea_surface_height_above_geoid','sea_surface_elevation',
    'sea_surface_height_above_reference_ellipsoid','sea_surface_height_above_sea_level',
    'sea_surface_height','water level']


    

In [5]:

    

print std_name_list
def name_in_list(cube):
    return cube.standard_name in std_name_list
constraint = iris.Constraint(cube_func=name_in_list)


    

In [6]:

    

a = iris.load_cube(url,constraint)


    

    




/home/local/python27_epd/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1216: UserWarning: Gracefully filling 'time' dimension coordinate masked points
  warnings.warn(msg.format(str(cf_coord_var.cf_name)))
/home/local/python27_epd/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1280: UserWarning: Failed to create 'time' dimension coordinate: The points array must be strictly monotonic.
Gracefully creating 'time' auxiliary coordinate instead.
  error=e_msg))

In [7]:

    

print a


    

    




water_surface_height_above_reference_datum / (ft) (-- : 50692)
     Auxiliary coordinates:
          time                                        x
     Scalar coordinates:
          latitude: 41.7207 degrees
          longitude: -73.939 degrees
     Attributes:
          Conventions: CF-1.6
          DODS.dimName: name_strlen
          DODS.strlen: 7
          _ChunkSize: 1
          comment: Datum is NAVD88
          date_created: 2014-03-16T05:00:00Z
          description: Water Elevation (NAVD88)
          featureType: timeSeries
          geospatial_lat_max: 41.7207
          geospatial_lat_min: 41.7207
          geospatial_lon_max: -73.939
          geospatial_lon_min: -73.939
          institution: HRECOS.
          keywords: Hudson, Weather, Wind, Hydro
          project: HRECOS
          publisher_email: devops@asascience.com
          publisher_name: RPS ASA on behalf of HRECOS.
          publisher_phone: (401) 789-6224
          publisher_url: http://www.asascience.com/
          references: http://www.HRECOS.org
          summary: Aggregated Dataset of HRECOS Station HRMARPH
          time_coverage_duration: P86400S
          time_coverage_end: 1394942400.0
          time_coverage_resolution: P360S
          time_coverage_start: 1394856000.0
          title: HRECOS Aggregated Station HRMARPH Data

In [8]:

    

a.convert_units('m')


    

In [9]:

    

a.coord('time')


    

    
Out[9]:





AuxCoord(array([  1.38483720e+09,   1.38483756e+09,   1.38483792e+09, ...,
         1.39502772e+09,   1.39502808e+09,   1.39502880e+09]), standard_name=u'time', units=Unit('seconds since 1970-01-01 00:00:00', calendar='gregorian'), long_name=u'time of measurement', attributes={'_CoordianteAxisType': 'Time', '_ChunkSize': 1})

In [10]:

    

def find_timevar(cube):
    """
    return the time variable from Iris. This is a workaround for
    Iris having problems with FMRC aggregations, which produce two time coordinates
    """
    try:
        cube.coord(axis='T').rename('time')
    except:
        pass
    timevar = cube.coord('time')
    return timevar


    

In [11]:

    

mod_name = a.attributes['title'][0:20]
r = shape(a)
timevar = find_timevar(a)
lat = a.coord(axis='Y').points
lon = a.coord(axis='X').points
jd = timevar.units.num2date(timevar.points)
istart = timevar.nearest_neighbour_index(timevar.units.date2num(jd_start))
istop = timevar.nearest_neighbour_index(timevar.units.date2num(jd_stop))


    

    




---------------------------------------------------------------------------
IllegalMonthError                         Traceback (most recent call last)
<ipython-input-11-8cf8e60dd282> in <module>()
      4 lat = a.coord(axis='Y').points
      5 lon = a.coord(axis='X').points
----> 6 jd = timevar.units.num2date(timevar.points)
      7 istart = timevar.nearest_neighbour_index(timevar.units.date2num(jd_start))
      8 istop = timevar.nearest_neighbour_index(timevar.units.date2num(jd_stop))

/home/local/python27_epd/lib/python2.7/site-packages/iris/unit.pyc in num2date(self, time_value)
   1972         """
   1973         cdf_utime = self.utime()
-> 1974         return cdf_utime.num2date(time_value)

/home/local/python27_epd/lib/python2.7/site-packages/netcdftime.pyc in num2date(self, time_value)
    772                             date.append(None)
    773                 else:
--> 774                     date = [DateFromJulianDay(j,self.calendar) for j in jd.flat]
    775             else:
    776                 if ismasked and mask.item():

/home/local/python27_epd/lib/python2.7/site-packages/netcdftime.pyc in DateFromJulianDay(JD, calendar)
    304     # if days exceeds number allowed in a month, flip to next month.
    305     # this fixes issue 75.
--> 306     daysinmonth = monthrange(year, month)[1]
    307     if days > daysinmonth:
    308         days = 1

/home/local/python27_epd/lib/python2.7/calendar.pyc in monthrange(year, month)
    118        year, month."""
    119     if not 1 <= month <= 12:
--> 120         raise IllegalMonthError(month)
    121     day1 = weekday(year, month, 1)
    122     ndays = mdays[month] + (month == February and isleap(year))

IllegalMonthError: bad month number -106757379390082; must be 1-12

In [12]:

    

t=timevar.points


    

In [13]:

    

t.max()


    

    
Out[13]:





9.969209968386869e+36

In [ ]: