1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import logging
7import os
8from collections import namedtuple
9import os.path as op
11import numpy as num
13from pyrocko import util, config
14from pyrocko import orthodrome as od
15from pyrocko.util import urlopen
17logger = logging.getLogger('pyrocko.dataset.geonames')
19GeoName = namedtuple('GeoName', '''
20geonameid
21name
22asciiname
23alt_names
24lat
25lon
26feature_class
27feature_code
28country_code
29alt_country_code
30admin1_code
31admin2_code
32admin3_code
33admin4_code
34population
35elevation
36dem
37timezone
38modification_date'''.split())
40GeoName2 = namedtuple('GeoName2', '''
41name
42asciiname
43lat
44lon
45population
46feature_code
47'''.split())
49base_url = 'https://mirror.pyrocko.org/download.geonames.org/export/dump/'
52def download_file(fn, dirpath):
53 url = base_url + '/' + fn
54 fpath = op.join(dirpath, fn)
55 logger.info('starting download of %s' % url)
57 util.ensuredirs(fpath)
58 f = urlopen(url)
59 fpath_tmp = fpath + '.%i.temp' % os.getpid()
60 g = open(fpath_tmp, 'wb')
61 while True:
62 data = f.read(1024)
63 if not data:
64 break
65 g.write(data)
67 g.close()
68 f.close()
70 os.rename(fpath_tmp, fpath)
72 logger.info('finished download of %s' % url)
75def positive_region(region):
76 west, east, south, north = [float(x) for x in region]
78 assert -180. - 360. <= west < 180.
79 assert -180. < east <= 180. + 360.
80 assert -90. <= south < 90.
81 assert -90. < north <= 90.
83 if east < west:
84 east += 360.
86 if west < -180.:
87 west += 360.
88 east += 360.
90 return (west, east, south, north)
93def ascii_str(u):
94 return u.encode('ascii', 'replace')
97def load(zfn, fn, minpop=1000000, region=None):
98 geonames_dir = config.config().geonames_dir
99 filepath = op.join(geonames_dir, zfn or fn)
100 if not os.path.exists(filepath):
101 download_file(zfn or fn, geonames_dir)
103 if region:
104 w, e, s, n = positive_region(region)
106 if zfn is not None:
107 import zipfile
108 z = zipfile.ZipFile(filepath, 'r')
109 f = z.open(fn, 'r')
110 else:
111 z = None
112 f = open(filepath, 'rb')
114 for line in f:
115 t = line.split(b'\t')
116 pop = int(t[14])
117 if minpop <= pop:
118 lat = float(t[4])
119 lon = float(t[5])
120 feature_code = str(t[7].decode('ascii'))
121 if not region or (
122 (w <= lon <= e or w <= lon + 360. <= e)
123 and (s <= lat <= n)):
125 yield GeoName2(
126 t[1].decode('utf8'),
127 str(t[2].decode('utf8').encode('ascii', 'replace')
128 .decode('ascii')),
129 lat, lon, pop, feature_code)
131 f.close()
132 if z is not None:
133 z.close()
136g_cities = {}
139def load_all_keep(zfn, fn, minpop=1000000, region=None, exclude=()):
140 if (zfn, fn) not in g_cities:
141 g_cities[zfn, fn] = list(load(zfn, fn, minpop=0))
143 if region:
144 w, e, s, n = positive_region(region)
145 return [c for c in g_cities[zfn, fn]
146 if (minpop <= c.population and
147 ((w <= c.lon <= e or w <= c.lon + 360. <= e)
148 and (s <= c.lat <= n)) and
149 c.feature_code not in exclude)]
151 else:
152 return [c for c in g_cities[zfn, fn] if (
153 minpop <= c.population and c.feature_code not in exclude)]
156def get_cities_region(region=None, minpop=0):
157 return load_all_keep(
158 'cities1000.zip', 'cities1000.txt',
159 region=region,
160 minpop=minpop,
161 exclude=('PPLX',))
164def get_cities_by_name(name):
165 cities = get_cities_region()
166 candidates = []
167 for c in cities:
168 if c.asciiname.lower() == name.lower():
169 candidates.append(c)
171 return candidates
174def get_cities(lat, lon, radius, minpop=0):
175 region = od.radius_to_region(lat, lon, radius)
176 cities = get_cities_region(region, minpop=minpop)
178 clats = num.array([c.lat for c in cities])
179 clons = num.array([c.lon for c in cities])
181 dists = od.distance_accurate50m_numpy(lat, lon, clats, clons)
182 order = num.argsort(dists)
183 cities_sorted = [cities[i] for i in order if dists[i] < radius]
184 return cities_sorted