|
# 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
|