1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

# http://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

# ---|P------/S----------~Lg---------- 

from __future__ import absolute_import, division 

 

try: 

from future.moves.urllib.request import urlopen 

except ImportError: 

from urllib2 import urlopen 

 

import time 

import re 

import logging 

import json 

 

from pyrocko import model, util 

from pyrocko.moment_tensor import MomentTensor, symmat6 

from .base_catalog import EarthquakeCatalog, NotFound 

 

logger = logging.getLogger('pyrocko.client.geofon') 

 

km = 1000. 

 

 

class Geofon(EarthquakeCatalog): 

'''Access the Geofon earthquake catalog ''' 

def __init__(self, get_moment_tensors=True): 

self.events = {} 

self._get_moment_tensors = get_moment_tensors 

 

def flush(self): 

self.events = {} 

 

def iter_event_names( 

self, 

time_range=None, 

nmax=1000, 

magmin=None, 

latmin=-90., 

latmax=90., 

lonmin=-180., 

lonmax=180.): 

 

logger.debug('In Geofon.iter_event_names(...)') 

 

dmin = time.strftime('%Y-%m-%d', time.gmtime(time_range[0])) 

dmax = time.strftime('%Y-%m-%d', time.gmtime(time_range[1]+24*60*60)) 

 

if magmin is None: 

magmin = '' 

else: 

magmin = '%g' % magmin 

 

ipage = 1 

while True: 

url = ('http://geofon.gfz-potsdam.de/eqinfo/list.php?' + '&'.join([ 

'page=%i' % ipage, 

'datemin=%s' % dmin, 

'datemax=%s' % dmax, 

'latmin=%g' % latmin, 

'latmax=%g' % latmax, 

'lonmin=%g' % lonmin, 

'lonmax=%g' % lonmax, 

'magmin=%s' % magmin, 

'fmt=geojson', 

'nmax=%i' % nmax])) 

 

logger.debug('Opening URL: %s' % url) 

page = urlopen(url).read() 

logger.debug('Received page (%i bytes)' % len(page)) 

events = self._parse_events_page(page) 

for ev in events: 

if ev.moment_tensor is True: 

ev.moment_tensor = self.get_mt(ev) 

 

if not events: 

break 

 

for ev in events: 

if time_range[0] <= ev.time and ev.time <= time_range[1]: 

self.events[ev.name] = ev 

yield ev.name 

 

ipage += 1 

 

def get_event(self, name): 

logger.debug('In Geofon.get_event("%s")' % name) 

 

if name not in self.events: 

url = 'http://geofon.gfz-potsdam.de/eqinfo/event.php' \ 

'?id=%s&fmt=geojson' % name 

logger.debug('Opening URL: %s' % url) 

page = urlopen(url).read() 

logger.debug('Received page (%i bytes)' % len(page)) 

 

try: 

ev = self._parse_event_page(page) 

self.events[name] = ev 

 

except NotFound: 

raise NotFound(url) # reraise with url 

 

ev = self.events[name] 

 

if ev.moment_tensor is True: 

ev.moment_tensor = self.get_mt(ev) 

 

return ev 

 

def get_mt(self, ev): 

syear = time.strftime('%Y', time.gmtime(ev.time)) 

url = 'http://geofon.gfz-potsdam.de/data/alerts/%s/%s/mt.txt' % ( 

syear, ev.name) 

logger.debug('Opening URL: %s' % url) 

page = urlopen(url).read() 

logger.debug('Received page (%i bytes)' % len(page)) 

 

return self._parse_mt_page(page) 

 

def _parse_mt_page(self, page): 

d = {} 

for k in 'Scale', 'Mrr', 'Mtt', 'Mpp', 'Mrt', 'Mrp', 'Mtp': 

r = k.encode('ascii')+br'\s*=?\s*(\S+)' 

m = re.search(r, page) 

if m: 

s = m.group(1).replace(b'10**', b'1e') 

d[k.lower()] = float(s) 

 

m = symmat6(*(d[x] for x in 'mrr mtt mpp mrt mrp mtp'.split())) 

m *= d['scale'] 

mt = MomentTensor(m_up_south_east=m) 

 

return mt 

 

def _parse_events_page(self, page, limit=None): 

j = json.loads(page.decode('utf-8')) 

events = [] 

for ifeature, feature in enumerate(j['features']): 

ev = self._json_feature_to_event(feature) 

events.append(ev) 

if limit and ifeature + 1 == limit: 

break 

 

return events 

 

def _parse_event_page(self, page): 

return self._parse_events_page(page, limit=1)[0] 

 

def _json_feature_to_event(self, feature): 

name = feature['id'] 

lon, lat, depth = feature['geometry']['coordinates'] 

depth *= 1000. 

properties = feature['properties'] 

magnitude = properties['mag'] 

magnitude_type = properties['magType'] 

region = properties['place'] 

tevent = util.str_to_time(properties['time'].replace('T', ' ')) 

 

if 'hasMT' in properties and properties['hasMT'] == 'yes' \ 

and self._get_moment_tensors: 

 

moment_tensor = True # flag for caller to query MT 

else: 

moment_tensor = None 

 

status = properties['status'][:1] 

tags = [] 

if status in 'AMC': 

tags.append('geofon_status:%s' % status) 

 

category = properties.get('evtype', '') 

if re.match(r'^[a-zA-Z0-9]+$', category): 

tags.append('geofon_category:%s' % category) 

 

ev = model.Event( 

lat=float(lat), 

lon=float(lon), 

time=tevent, 

name=name, 

depth=float(depth), 

magnitude=float(magnitude), 

magnitude_type=str(magnitude_type), 

region=str(region), 

moment_tensor=moment_tensor, 

catalog='GEOFON', 

tags=tags) 

 

return ev