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