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