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