1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3"""
4Provides the Station class.
5
6:copyright:
7    Lion Krischer (krischer@geophysik.uni-muenchen.de), 2013
8:license:
9    GNU Lesser General Public License, Version 3
10    (https://www.gnu.org/copyleft/lesser.html)
11"""
12from __future__ import (absolute_import, division, print_function,
13                        unicode_literals)
14from future.builtins import *  # NOQA
15from future.utils import python_2_unicode_compatible
16
17import copy
18import fnmatch
19import warnings
20
21import numpy as np
22
23from obspy import UTCDateTime
24from obspy.core.util.obspy_types import (ObsPyException, ZeroSamplingRate,
25                                         FloatWithUncertaintiesAndUnit)
26from obspy.geodetics import inside_geobounds
27
28from .util import (BaseNode, Equipment, Operator, Distance, Latitude,
29                   Longitude, _unified_content_strings, _textwrap, Site)
30
31
32@python_2_unicode_compatible
33class Station(BaseNode):
34    """
35    From the StationXML definition:
36        This type represents a Station epoch. It is common to only have a
37        single station epoch with the station's creation and termination dates
38        as the epoch start and end dates.
39    """
40    def __init__(self, code, latitude, longitude, elevation, channels=None,
41                 site=None, vault=None, geology=None, equipments=None,
42                 operators=None, creation_date=None, termination_date=None,
43                 total_number_of_channels=None,
44                 selected_number_of_channels=None, description=None,
45                 comments=None, start_date=None, end_date=None,
46                 restricted_status=None, alternate_code=None,
47                 historical_code=None, data_availability=None,
48                 identifiers=None, water_level=None, source_id=None):
49        """
50        :type channels: list of :class:`~obspy.core.inventory.channel.Channel`
51        :param channels: All channels belonging to this station.
52        :type site: :class:`~obspy.core.inventory.util.Site`
53        :param site: The lexical description of the site
54        :type latitude: :class:`~obspy.core.inventory.util.Latitude`
55        :param latitude: The latitude of the station
56        :type longitude: :class:`~obspy.core.inventory.util.Longitude`
57        :param longitude: The longitude of the station
58        :param elevation: The elevation of the station in meter.
59        :param site: These fields describe the location of the station using
60            geopolitical entities (country, city, etc.).
61        :param vault: Type of vault, e.g. WWSSN, tunnel, transportable array,
62            etc
63        :param geology: Type of rock and/or geologic formation.
64        :type equipments: list of :class:`~obspy.core.inventory.util.Equipment`
65        :param equipments: Equipment used by all channels at a station.
66        :type operators: list of :class:`~obspy.core.inventory.util.Operator`
67        :param operators: An operating agency and associated contact persons.
68        :type creation_date: :class:`~obspy.core.utcdatetime.UTCDateTime`
69        :param creation_date: Date and time (UTC) when the station was first
70            installed
71        :type termination_date: :class:`~obspy.core.utcdatetime.UTCDateTime`,
72            optional
73        :param termination_date: Date and time (UTC) when the station was
74            terminated or will be terminated. A blank value should be assumed
75            to mean that the station is still active.
76        :type total_number_of_channels: int
77        :param total_number_of_channels: Total number of channels recorded at
78            this station.
79        :type selected_number_of_channels: int
80        :param selected_number_of_channels: Number of channels recorded at this
81            station and selected by the query that produced this document.
82        :type external_references: list of
83            :class:`~obspy.core.inventory.util.ExternalReference`
84        :param external_references: URI of any type of external report, such as
85            IRIS data reports or dataless SEED volumes.
86        :type description: str
87        :param description: A description of the resource
88        :type comments: list of :class:`~obspy.core.inventory.util.Comment`
89        :param comments: An arbitrary number of comments to the resource
90        :type start_date: :class:`~obspy.core.utcdatetime.UTCDateTime`,
91            optional
92        :param start_date: The start date of the resource
93        :type end_date: :class:`~obspy.core.utcdatetime.UTCDateTime`
94        :param end_date: The end date of the resource
95        :type restricted_status: str
96        :param restricted_status: The restriction status
97        :type alternate_code: str
98        :param alternate_code: A code used for display or association,
99            alternate to the SEED-compliant code.
100        :type historical_code: str
101        :param historical_code: A previously used code if different from the
102            current code.
103        :type data_availability: :class:`~obspy.station.util.DataAvailability`
104        :param data_availability: Information about time series availability
105            for the station.
106        :type identifiers: list of str, optional
107        :param identifiers: Persistent identifiers for network/station/channel
108            (schema version >=1.1). URIs are in general composed of a 'scheme'
109            and a 'path' (optionally with additional components), the two of
110            which separated by a colon.
111        :type water_level: float, optional
112        :param water_level: Elevation of the water surface in meters for
113            underwater sites, where 0 is sea level. (schema version >=1.1)
114        :type source_id: str, optional
115        :param source_id: A data source identifier in URI form
116            (schema version >=1.1). URIs are in general composed of a 'scheme'
117            and a 'path' (optionally with additional components), the two of
118            which separated by a colon.
119        """
120        self.latitude = latitude
121        self.longitude = longitude
122        self.elevation = elevation
123        self.channels = channels or []
124        self.site = site if site is not None else Site()
125        self.vault = vault
126        self.geology = geology
127        self.equipments = equipments or []
128        self.operators = operators or []
129        self.creation_date = creation_date
130        self.termination_date = termination_date
131        self.total_number_of_channels = total_number_of_channels
132        self.selected_number_of_channels = selected_number_of_channels
133        self.external_references = []
134        self.water_level = water_level
135        super(Station, self).__init__(
136            code=code, description=description, comments=comments,
137            start_date=start_date, end_date=end_date,
138            restricted_status=restricted_status, alternate_code=alternate_code,
139            historical_code=historical_code,
140            data_availability=data_availability, identifiers=identifiers,
141            source_id=source_id)
142
143    @property
144    def total_number_of_channels(self):
145        return self._total_number_of_channels
146
147    @total_number_of_channels.setter
148    def total_number_of_channels(self, value):
149        if value is not None and value < 0:
150            msg = "total_number_of_channels cannot be negative."
151            raise ValueError(msg)
152        self._total_number_of_channels = value
153
154    @property
155    def selected_number_of_channels(self):
156        return self._selected_number_of_channels
157
158    @selected_number_of_channels.setter
159    def selected_number_of_channels(self, value):
160        if value is not None and value < 0:
161            msg = "selected_number_of_channels cannot be negative."
162            raise ValueError(msg)
163        self._selected_number_of_channels = value
164
165    def __str__(self):
166        contents = self.get_contents()
167        ret = ("Station {station_name}\n"
168               "\tStation Code: {station_code}\n"
169               "\tChannel Count: {selected}/{total} (Selected/Total)\n"
170               "\t{start_date} - {end_date}\n"
171               "\tAccess: {restricted} {alternate_code}{historical_code}\n"
172               "\tLatitude: {lat:.2f}, Longitude: {lng:.2f}, "
173               "Elevation: {elevation:.1f} m\n")
174        ret = ret.format(
175            station_name=contents["stations"][0],
176            station_code=self.code,
177            selected=self.selected_number_of_channels,
178            total=self.total_number_of_channels,
179            start_date=str(self.start_date),
180            end_date=str(self.end_date) if self.end_date else "",
181            restricted=self.restricted_status,
182            lat=self.latitude, lng=self.longitude, elevation=self.elevation,
183            alternate_code="Alternate Code: %s " % self.alternate_code if
184            self.alternate_code else "",
185            historical_code="historical Code: %s " % self.historical_code if
186            self.historical_code else "")
187        ret += "\tAvailable Channels:\n"
188        ret += "\n".join(_textwrap(
189            ", ".join(_unified_content_strings(contents["channels"])),
190            initial_indent="\t\t", subsequent_indent="\t\t",
191            expand_tabs=False))
192        return ret
193
194    def _repr_pretty_(self, p, cycle):
195        p.text(str(self))
196
197    def __getitem__(self, index):
198        return self.channels[index]
199
200    def __len__(self):
201        return len(self.channels)
202
203    def get_contents(self):
204        """
205        Returns a dictionary containing the contents of the object.
206
207        .. rubric:: Example
208
209        >>> from obspy import read_inventory
210        >>> example_filename = "/path/to/IRIS_single_channel_with_response.xml"
211        >>> inventory = read_inventory(example_filename)
212        >>> station = inventory.networks[0].stations[0]
213        >>> station.get_contents()  # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
214        {...}
215        >>> for (k, v) in sorted(station.get_contents().items()):
216        ...     print(k, v[0])
217        channels ANMO.10.BHZ
218        stations ANMO (Albuquerque, New Mexico, USA)
219        """
220        site_name = None
221        if self.site and self.site.name:
222            site_name = self.site.name
223        desc = "%s%s" % (self.code, " (%s)" % (site_name if site_name else ""))
224        content_dict = {"stations": [desc], "channels": []}
225
226        for channel in self.channels:
227            content_dict["channels"].append(
228                "%s.%s.%s" % (self.code, channel.location_code, channel.code))
229        return content_dict
230
231    @property
232    def operators(self):
233        return self._operators
234
235    @operators.setter
236    def operators(self, value):
237        if not hasattr(value, "__iter__"):
238            msg = "Operators needs to be an iterable, e.g. a list."
239            raise ValueError(msg)
240        # make sure to unwind actual iterators, or the just might get exhausted
241        # at some point
242        operators = [operator for operator in value]
243        if any([not isinstance(x, Operator) for x in operators]):
244            msg = "Operators can only contain Operator objects."
245            raise ValueError(msg)
246        self._operators = operators
247
248    @property
249    def equipments(self):
250        return self._equipments
251
252    @equipments.setter
253    def equipments(self, value):
254        if not hasattr(value, "__iter__"):
255            msg = "equipments needs to be an iterable, e.g. a list."
256            raise ValueError(msg)
257        # make sure to unwind actual iterators, or the just might get exhausted
258        # at some point
259        equipments = [equipment for equipment in value]
260        if any([not isinstance(x, Equipment) for x in equipments]):
261            msg = "equipments can only contain Equipment objects."
262            raise ValueError(msg)
263        self._equipments = equipments
264        # if value is None or isinstance(value, Equipment):
265        #    self._equipment = value
266        # elif isinstance(value, dict):
267        #    self._equipment = Equipment(**value)
268        # else:
269        #    msg = ("equipment needs to be be of type
270        # obspy.core.inventory.Equipment "
271        #        "or contain a dictionary with values suitable for "
272        #        "initialization.")
273        #    raise ValueError(msg)
274
275    @property
276    def creation_date(self):
277        return self._creation_date
278
279    @creation_date.setter
280    def creation_date(self, value):
281        if value is None:
282            self._creation_date = None
283            return
284        if not isinstance(value, UTCDateTime):
285            value = UTCDateTime(value)
286        self._creation_date = value
287
288    @property
289    def termination_date(self):
290        return self._termination_date
291
292    @termination_date.setter
293    def termination_date(self, value):
294        if value is not None and not isinstance(value, UTCDateTime):
295            value = UTCDateTime(value)
296        self._termination_date = value
297
298    @property
299    def external_references(self):
300        return self._external_references
301
302    @external_references.setter
303    def external_references(self, value):
304        if not hasattr(value, "__iter__"):
305            msg = "external_references needs to be iterable, e.g. a list."
306            raise ValueError(msg)
307        self._external_references = value
308
309    @property
310    def longitude(self):
311        return self._longitude
312
313    @longitude.setter
314    def longitude(self, value):
315        if isinstance(value, Longitude):
316            self._longitude = value
317        else:
318            self._longitude = Longitude(value)
319
320    @property
321    def latitude(self):
322        return self._latitude
323
324    @latitude.setter
325    def latitude(self, value):
326        if isinstance(value, Latitude):
327            self._latitude = value
328        else:
329            self._latitude = Latitude(value)
330
331    @property
332    def elevation(self):
333        return self._elevation
334
335    @elevation.setter
336    def elevation(self, value):
337        if isinstance(value, Distance):
338            self._elevation = value
339        else:
340            self._elevation = Distance(value)
341
342    @property
343    def water_level(self):
344        return self._water_level
345
346    @water_level.setter
347    def water_level(self, value):
348        if value is None:
349            self._water_level = None
350        elif isinstance(value, FloatWithUncertaintiesAndUnit):
351            self._water_level = value
352        else:
353            self._water_level = FloatWithUncertaintiesAndUnit(value)
354
355    def select(self, location=None, channel=None, time=None, starttime=None,
356               endtime=None, sampling_rate=None, minlatitude=None,
357               maxlatitude=None, minlongitude=None, maxlongitude=None,
358               latitude=None, longitude=None, minradius=None, maxradius=None):
359        r"""
360        Returns the :class:`Station` object with only the
361        :class:`~obspy.core.inventory.channel.Channel`\ s that match the given
362        criteria (e.g. all channels with ``channel="EHZ"``).
363
364        .. warning::
365            The returned object is based on a shallow copy of the original
366            object. That means that modifying any mutable child elements will
367            also modify the original object
368            (see https://docs.python.org/2/library/copy.html).
369            Use :meth:`copy()` afterwards to make a new copy of the data in
370            memory.
371
372        .. rubric:: Example
373
374        >>> from obspy import read_inventory, UTCDateTime
375        >>> sta = read_inventory()[0][0]
376        >>> t = UTCDateTime(2008, 7, 1, 12)
377        >>> sta = sta.select(channel="[LB]HZ", time=t)
378        >>> print(sta)  # doctest: +NORMALIZE_WHITESPACE
379        Station FUR (Fuerstenfeldbruck, Bavaria, GR-Net)
380            Station Code: FUR
381            Channel Count: None/None (Selected/Total)
382            2006-12-16T00:00:00.000000Z -
383            Access: None
384            Latitude: 48.16, Longitude: 11.28, Elevation: 565.0 m
385            Available Channels:
386                FUR..BHZ, FUR..LHZ
387
388        The `location` and `channel` selection criteria  may also contain UNIX
389        style wildcards (e.g. ``*``, ``?``, ...; see
390        :func:`~fnmatch.fnmatch`).
391
392        :type location: str
393        :param location: Potentially wildcarded location code. If not given,
394            all location codes will be accepted.
395        :type channel: str
396        :param channel: Potentially wildcarded channel code. If not given,
397            all channel codes will be accepted.
398        :type time: :class:`~obspy.core.utcdatetime.UTCDateTime`
399        :param time: Only include channels active at given point in time.
400        :type starttime: :class:`~obspy.core.utcdatetime.UTCDateTime`
401        :param starttime: Only include channels active at or after given point
402            in time (i.e. channels ending before given time will not be shown).
403        :type endtime: :class:`~obspy.core.utcdatetime.UTCDateTime`
404        :param endtime: Only include channels active before or at given point
405            in time (i.e. channels starting after given time will not be
406            shown).
407        :type sampling_rate: float
408        :param sampling_rate: Only include channels whose sampling rate
409            matches the given sampling rate, in Hz (within absolute tolerance
410            of 1E-8 Hz and relative tolerance of 1E-5)
411        :type minlatitude: float
412        :param minlatitude: Only include channels with a latitude larger than
413            the specified minimum.
414        :type maxlatitude: float
415        :param maxlatitude: Only include channels with a latitude smaller than
416            the specified maximum.
417        :type minlongitude: float
418        :param minlongitude: Only include channels with a longitude larger than
419            the specified minimum.
420        :type maxlongitude: float
421        :param maxlongitude: Only include channels with a longitude smaller
422            than the specified maximum.
423        :type latitude: float
424        :param latitude: Specify the latitude to be used for a radius
425            selection.
426        :type longitude: float
427        :param longitude: Specify the longitude to be used for a radius
428            selection.
429        :type minradius: float
430        :param minradius: Only include channels within the specified
431            minimum number of degrees from the geographic point defined by the
432            latitude and longitude parameters.
433        :type maxradius: float
434        :param maxradius: Only include channels within the specified
435            maximum number of degrees from the geographic point defined by the
436            latitude and longitude parameters.
437        """
438        channels = []
439        for cha in self.channels:
440            # skip if any given criterion is not matched
441            if location is not None:
442                if not fnmatch.fnmatch(cha.location_code.upper(),
443                                       location.upper()):
444                    continue
445            if channel is not None:
446                if not fnmatch.fnmatch(cha.code.upper(),
447                                       channel.upper()):
448                    continue
449            if sampling_rate is not None:
450                if cha.sample_rate is None:
451                    msg = ("Omitting channel that has no sampling rate "
452                           "specified.")
453                    warnings.warn(msg)
454                    continue
455                if not np.allclose(float(sampling_rate), cha.sample_rate,
456                                   rtol=1E-5, atol=1E-8):
457                    continue
458            if any([t is not None for t in (time, starttime, endtime)]):
459                if not cha.is_active(time=time, starttime=starttime,
460                                     endtime=endtime):
461                    continue
462            geo_filters = dict(
463                minlatitude=minlatitude, maxlatitude=maxlatitude,
464                minlongitude=minlongitude, maxlongitude=maxlongitude,
465                latitude=latitude, longitude=longitude, minradius=minradius,
466                maxradius=maxradius)
467            if any(value is not None for value in geo_filters.values()):
468                if not inside_geobounds(cha, **geo_filters):
469                    continue
470
471            channels.append(cha)
472        sta = copy.copy(self)
473        sta.channels = channels
474        return sta
475
476    def plot(self, min_freq, output="VEL", location="*", channel="*",
477             time=None, starttime=None, endtime=None, axes=None,
478             unwrap_phase=False, plot_degrees=False, show=True, outfile=None):
479        """
480        Show bode plot of instrument response of all (or a subset of) the
481        station's channels.
482
483        :type min_freq: float
484        :param min_freq: Lowest frequency to plot.
485        :type output: str
486        :param output: Output units. One of:
487
488                ``"DISP"``
489                    displacement, output unit is meters
490                ``"VEL"``
491                    velocity, output unit is meters/second
492                ``"ACC"``
493                    acceleration, output unit is meters/second**2
494
495        :type location: str
496        :param location: Only plot matching channels. Accepts UNIX style
497            patterns and wildcards (e.g. ``"BH*"``, ``"BH?"``, ``"*Z"``,
498            ``"[LB]HZ"``; see :func:`~fnmatch.fnmatch`)
499        :type channel: str
500        :param channel: Only plot matching channels. Accepts UNIX style
501            patterns and wildcards (e.g. ``"BH*"``, ``"BH?"``, ``"*Z"``,
502            ``"[LB]HZ"``; see :func:`~fnmatch.fnmatch`)
503        :param time: Only show channels active at given point in time.
504        :type starttime: :class:`~obspy.core.utcdatetime.UTCDateTime`
505        :param starttime: Only show channels active at or after given point in
506            time (i.e. channels ending before given time will not be shown).
507        :type endtime: :class:`~obspy.core.utcdatetime.UTCDateTime`
508        :param endtime: Only show channels active before or at given point in
509            time (i.e. channels starting after given time will not be shown).
510        :type axes: list of 2 :class:`matplotlib.axes.Axes`
511        :param axes: List/tuple of two axes instances to plot the
512            amplitude/phase spectrum into. If not specified, a new figure is
513            opened.
514        :type unwrap_phase: bool
515        :param unwrap_phase: Set optional phase unwrapping using NumPy.
516        :type plot_degrees: bool
517        :param plot_degrees: if ``True`` plot bode in degrees
518        :type show: bool
519        :param show: Whether to show the figure after plotting or not. Can be
520            used to do further customization of the plot before showing it.
521        :type outfile: str
522        :param outfile: Output file path to directly save the resulting image
523            (e.g. ``"/tmp/image.png"``). Overrides the ``show`` option, image
524            will not be displayed interactively. The given path/file name is
525            also used to automatically determine the output format. Supported
526            file formats depend on your matplotlib backend.  Most backends
527            support png, pdf, ps, eps and svg. Defaults to ``None``.
528
529        .. rubric:: Basic Usage
530
531        >>> from obspy import read_inventory
532        >>> sta = read_inventory()[0][0]
533        >>> sta.plot(0.001, output="VEL", channel="*Z")  # doctest: +SKIP
534
535        .. plot::
536
537            from obspy import read_inventory
538            sta = read_inventory()[0][0]
539            sta.plot(0.001, output="VEL", channel="*Z")
540        """
541        import matplotlib.pyplot as plt
542
543        if axes is not None:
544            ax1, ax2 = axes
545            fig = ax1.figure
546        else:
547            fig = plt.figure()
548            ax1 = fig.add_subplot(211)
549            ax2 = fig.add_subplot(212, sharex=ax1)
550
551        matching = self.select(location=location, channel=channel, time=time,
552                               starttime=starttime, endtime=endtime)
553
554        for cha in matching.channels:
555            try:
556                cha.plot(min_freq=min_freq, output=output, axes=(ax1, ax2),
557                         label=".".join((self.code, cha.location_code,
558                                         cha.code)),
559                         unwrap_phase=unwrap_phase, plot_degrees=plot_degrees,
560                         show=False, outfile=None)
561            except ZeroSamplingRate:
562                msg = ("Skipping plot of channel with zero "
563                       "sampling rate:\n%s")
564                warnings.warn(msg % str(cha), UserWarning)
565            except ObsPyException as e:
566                msg = "Skipping plot of channel (%s):\n%s"
567                warnings.warn(msg % (str(e), str(cha)), UserWarning)
568
569        # final adjustments to plot if we created the figure in here
570        if not axes:
571            from obspy.core.inventory.response import _adjust_bode_plot_figure
572            _adjust_bode_plot_figure(fig, plot_degrees=plot_degrees,
573                                     show=False)
574
575        if outfile:
576            fig.savefig(outfile)
577        else:
578            if show:
579                plt.show()
580
581        return fig
582
583
584if __name__ == '__main__':
585    import doctest
586    doctest.testmod(exclude_empty=True)
587