1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import math
8import logging
9import os.path as op
11import numpy as num
13from . import tile
14from ..util import get_download_callback
15from pyrocko import util
17try:
18 range = xrange
19except NameError:
20 pass
22logger = logging.getLogger('pyrocko.dataset.topo.dataset')
25class TiledGlobalDataset(object):
27 def __init__(self, name, nx, ny, ntx, nty, dtype, data_dir=None,
28 citation=None, region=None):
30 # dataset geometry (including overlap/endpoints)
31 self.nx = int(nx)
32 self.ny = int(ny)
33 self.xmin = -180.
34 self.xmax = 180.
35 self.ymin = -90.
36 self.ymax = 90.
37 self.dx = (self.xmax-self.xmin) / (self.nx - 1)
38 self.dy = (self.ymax-self.ymin) / (self.ny - 1)
40 # tile geometry (including overlap)
41 self.ntx = int(ntx)
42 self.nty = int(nty)
43 self.stx = (self.ntx - 1) * self.dx
44 self.sty = (self.nty - 1) * self.dy
46 self.ntilesx = (self.nx - 1) // (self.ntx - 1)
47 self.ntilesy = (self.ny - 1) // (self.nty - 1)
49 self.name = name
50 self.dtype = dtype
51 self.data_dir = data_dir
52 self.citation = citation
53 if region is not None:
54 self.region = self.positive_region(region)
55 else:
56 self.region = None
58 def covers(self, region):
59 if self.region is None:
60 return True
62 a = self.region
63 b = self.positive_region(region)
64 return ((
65 (a[0] <= b[0] and b[1] <= b[1]) or
66 (a[0] <= b[0]+360. and b[1]+360 <= a[1]))
67 and a[2] <= b[2] and b[3] <= a[3])
69 def is_suitable(self, region, dmin, dmax):
70 d = 360. / (self.nx - 1)
71 return self.covers(region) and dmin <= d <= dmax
73 def download_file(self, url, fpath, username=None, password=None):
74 # TODO: Add logstring
75 util.download_file(
76 url, fpath, username, password,
77 status_callback=get_download_callback(
78 'Downloading %s topography...' % self.__class__.__name__))
80 def x(self):
81 return self.xmin + num.arange(self.nx) * self.dx
83 def y(self):
84 return self.ymin + num.arange(self.ny) * self.dy
86 def positive_region(self, region):
87 xmin, xmax, ymin, ymax = [float(x) for x in region]
89 assert -180. - 360. <= xmin < 180.
90 assert -180. < xmax <= 180. + 360.
91 assert -90. <= ymin < 90.
92 assert -90. < ymax <= 90.
94 if xmax < xmin:
95 xmax += 360.
97 if xmin < -180.:
98 xmin += 360.
99 xmax += 360.
101 return (xmin, xmax, ymin, ymax)
103 def tile_indices(self, region):
104 xmin, xmax, ymin, ymax = self.positive_region(region)
105 itxmin = int(math.floor((xmin - self.xmin) / self.stx))
106 itxmax = int(math.ceil((xmax - self.xmin) / self.stx))
107 if itxmin == itxmax:
108 itxmax += 1
109 itymin = int(math.floor((ymin - self.ymin) / self.sty))
110 itymax = int(math.ceil((ymax - self.ymin) / self.sty))
111 if itymin == itymax:
112 if itymax != self.ntilesy:
113 itymax += 1
114 else:
115 itymin -= 1
117 indices = []
118 for ity in range(itymin, itymax):
119 for itx in range(itxmin, itxmax):
120 indices.append((itx % self.ntilesx, ity))
122 return indices
124 def get_tile(self, itx, ity):
125 return None
127 def get(self, region):
128 if len(region) == 2:
129 x, y = region
130 t = self.get((x, x, y, y))
131 if t is not None:
132 return t.get(x, y)
133 else:
134 return None
136 indices = self.tile_indices(region)
137 tiles = []
138 for itx, ity in indices:
139 t = self.get_tile(itx, ity)
140 if t:
141 tiles.append(t)
143 return tile.combine(tiles, region)
145 def get_with_repeat(self, region):
146 xmin, xmax, ymin, ymax = region
147 ymin2 = max(-90., ymin)
148 ymax2 = min(90., ymax)
149 region2 = xmin, xmax, ymin2, ymax2
150 t = self.get(region2)
151 if t is not None and region2 != region:
152 t.yextend_with_repeat(ymin, ymax)
154 return t
157class DecimatedTiledGlobalDataset(TiledGlobalDataset):
159 def __init__(self, name, base, ndeci, data_dir=None, ntx=None, nty=None):
161 assert ndeci % 2 == 0
162 assert (base.nx - 1) % ndeci == 0
163 assert (base.ny - 1) % ndeci == 0
165 nx = (base.nx - 1) // ndeci + 1
166 ny = (base.ny - 1) // ndeci + 1
168 if ntx is None:
169 ntx = base.ntx
171 if nty is None:
172 nty = base.nty
174 assert (nx - 1) % (ntx - 1) == 0
175 assert (ny - 1) % (nty - 1) == 0
177 if data_dir is None:
178 data_dir = op.join(base.data_dir, 'decimated_%i' % ndeci)
180 TiledGlobalDataset.__init__(self, name, nx, ny, ntx, nty, base.dtype,
181 data_dir=data_dir, citation=base.citation,
182 region=base.region)
184 self.base = base
185 self.ndeci = ndeci
187 def make_tile(self, itx, ity, fpath):
188 nh = self.ndeci // 2
189 xmin = self.xmin + itx*self.stx - self.base.dx * nh
190 xmax = self.xmin + (itx+1)*self.stx + self.base.dx * nh
191 ymin = self.ymin + ity*self.sty - self.base.dy * nh
192 ymax = self.ymin + (ity+1)*self.sty + self.base.dy * nh
194 t = self.base.get_with_repeat((xmin, xmax, ymin, ymax))
195 if t is not None:
196 t.decimate(self.ndeci)
198 util.ensuredirs(fpath)
199 with open(fpath, 'w') as f:
200 if t is not None:
201 t.data.tofile(f)
203 def make_if_needed(self, itx, ity):
204 assert 0 <= itx < self.ntilesx
205 assert 0 <= ity < self.ntilesy
207 fn = '%02i.%02i.bin' % (ity, itx)
208 fpath = op.join(self.data_dir, fn)
209 if not op.exists(fpath):
210 logger.info('making decimated tile: %s (%s)' % (fn, self.name))
211 self.make_tile(itx, ity, fpath)
213 def get_tile(self, itx, ity):
214 assert 0 <= itx < self.ntilesx
215 assert 0 <= ity < self.ntilesy
217 self.make_if_needed(itx, ity)
219 fn = '%02i.%02i.bin' % (ity, itx)
220 fpath = op.join(self.data_dir, fn)
221 with open(fpath, 'r') as f:
222 data = num.fromfile(f, dtype=self.dtype)
224 if data.size == 0:
225 return None
227 assert data.size == self.ntx*self.nty
228 data = data.reshape((self.ntx, self.nty))
230 return tile.Tile(
231 self.xmin + itx*self.stx,
232 self.ymin + ity*self.sty,
233 self.dx, self.dx, data)
235 def make_all_missing(self):
236 for ity in range(self.ntilesy):
237 for itx in range(self.ntilesx):
238 self.make_if_needed(itx, ity)