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

# http://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

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

from __future__ import absolute_import 

from future import standard_library 

standard_library.install_aliases() # noqa 

 

import logging 

import os 

from collections import namedtuple 

import os.path as op 

 

import numpy as num 

 

from pyrocko import util, config 

from pyrocko import orthodrome as od 

 

logger = logging.getLogger('pyrocko.dataset.geonames') 

 

GeoName = namedtuple('GeoName', ''' 

geonameid 

name 

asciiname 

alt_names 

lat 

lon 

feature_class 

feature_code 

country_code 

alt_country_code 

admin1_code 

admin2_code 

admin3_code 

admin4_code 

population 

elevation 

dem 

timezone 

modification_date'''.split()) 

 

GeoName2 = namedtuple('GeoName2', ''' 

name 

asciiname 

lat 

lon 

population 

feature_code 

'''.split()) 

 

base_url = 'https://mirror.pyrocko.org/download.geonames.org/export/dump/' 

 

 

def download_file(fn, dirpath): 

import urllib.request 

url = base_url + '/' + fn 

fpath = op.join(dirpath, fn) 

logger.info('starting download of %s' % url) 

 

util.ensuredirs(fpath) 

f = urllib.request.urlopen(url) 

fpath_tmp = fpath + '.%i.temp' % os.getpid() 

g = open(fpath_tmp, 'wb') 

while True: 

data = f.read(1024) 

if not data: 

break 

g.write(data) 

 

g.close() 

f.close() 

 

os.rename(fpath_tmp, fpath) 

 

logger.info('finished download of %s' % url) 

 

 

def positive_region(region): 

west, east, south, north = [float(x) for x in region] 

 

assert -180. - 360. <= west < 180. 

assert -180. < east <= 180. + 360. 

assert -90. <= south < 90. 

assert -90. < north <= 90. 

 

if east < west: 

east += 360. 

 

if west < -180.: 

west += 360. 

east += 360. 

 

return (west, east, south, north) 

 

 

def ascii_str(u): 

return u.encode('ascii', 'replace') 

 

 

def load(zfn, fn, minpop=1000000, region=None): 

geonames_dir = config.config().geonames_dir 

filepath = op.join(geonames_dir, zfn or fn) 

if not os.path.exists(filepath): 

download_file(zfn or fn, geonames_dir) 

 

if region: 

w, e, s, n = positive_region(region) 

 

if zfn is not None: 

import zipfile 

z = zipfile.ZipFile(filepath, 'r') 

f = z.open(fn, 'r') 

else: 

z = None 

f = open(filepath, 'rb') 

 

for line in f: 

t = line.split(b'\t') 

pop = int(t[14]) 

if minpop <= pop: 

lat = float(t[4]) 

lon = float(t[5]) 

feature_code = str(t[7].decode('ascii')) 

if not region or ( 

(w <= lon <= e or w <= lon + 360. <= e) 

and (s <= lat <= n)): 

 

yield GeoName2( 

t[1].decode('utf8'), 

str(t[2].decode('utf8').encode('ascii', 'replace') 

.decode('ascii')), 

lat, lon, pop, feature_code) 

 

f.close() 

if z is not None: 

z.close() 

 

 

g_cities = {} 

 

 

def load_all_keep(zfn, fn, minpop=1000000, region=None, exclude=()): 

if (zfn, fn) not in g_cities: 

g_cities[zfn, fn] = list(load(zfn, fn, minpop=0)) 

 

if region: 

w, e, s, n = positive_region(region) 

return [c for c in g_cities[zfn, fn] 

if (minpop <= c.population and 

((w <= c.lon <= e or w <= c.lon + 360. <= e) 

and (s <= c.lat <= n)) and 

c.feature_code not in exclude)] 

 

else: 

return [c for c in g_cities[zfn, fn] if ( 

minpop <= c.population and c.feature_code not in exclude)] 

 

 

def get_cities_region(region=None, minpop=0): 

return load_all_keep( 

'cities1000.zip', 'cities1000.txt', 

region=region, 

minpop=minpop, 

exclude=('PPLX',)) 

 

 

def get_cities_by_name(name): 

cities = get_cities_region() 

candidates = [] 

for c in cities: 

if c.asciiname.lower() == name.lower(): 

candidates.append(c) 

 

return candidates 

 

 

def get_cities(lat, lon, radius, minpop=0): 

region = od.radius_to_region(lat, lon, radius) 

cities = get_cities_region(region, minpop=minpop) 

 

clats = num.array([c.lat for c in cities]) 

clons = num.array([c.lon for c in cities]) 

 

dists = od.distance_accurate50m_numpy(lat, lon, clats, clons) 

order = num.argsort(dists) 

cities_sorted = [cities[i] for i in order if dists[i] < radius] 

return cities_sorted