1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import logging
8import os.path as op
10import numpy as num
12from . import tile
13from ..util import get_download_callback
14from pyrocko import util
16try:
17 range = xrange
18except NameError:
19 pass
21logger = logging.getLogger('pyrocko.dataset.topo.dataset')
24class TiledGlobalDataset(object):
26 def __init__(self, name, nx, ny, ntx, nty, dtype, data_dir=None,
27 citation=None, region=None):
29 # dataset geometry (including overlap/endpoints)
30 self.nx = int(nx)
31 self.ny = int(ny)
32 self.xmin = -180.
33 self.xmax = 180.
34 self.ymin = -90.
35 self.ymax = 90.
36 self.dx = (self.xmax-self.xmin) / (self.nx - 1)
37 self.dy = (self.ymax-self.ymin) / (self.ny - 1)
39 # tile geometry (including overlap)
40 self.ntx = int(ntx)
41 self.nty = int(nty)
42 self.stx = (self.ntx - 1) * self.dx
43 self.sty = (self.nty - 1) * self.dy
45 self.ntilesx = (self.nx - 1) // (self.ntx - 1)
46 self.ntilesy = (self.ny - 1) // (self.nty - 1)
48 self.name = name
49 self.dtype = dtype
50 self.data_dir = data_dir
51 self.citation = citation
52 if region is not None:
53 self.region = self.positive_region(region)
54 else:
55 self.region = None
57 def covers(self, region):
58 if self.region is None:
59 return True
61 a = self.region
62 b = self.positive_region(region)
63 return ((
64 (a[0] <= b[0] and b[1] <= b[1]) or
65 (a[0] <= b[0]+360. and b[1]+360 <= a[1]))
66 and a[2] <= b[2] and b[3] <= a[3])
68 def is_suitable(self, region, dmin, dmax):
69 d = 360. / (self.nx - 1)
70 return self.covers(region) and dmin <= d <= dmax
72 def download_file(self, url, fpath, username=None, password=None):
73 # TODO: Add logstring
74 util.download_file(
75 url, fpath, username, password,
76 status_callback=get_download_callback(
77 'Downloading %s topography...' % self.__class__.__name__))
79 def x(self):
80 return self.xmin + num.arange(self.nx) * self.dx
82 def y(self):
83 return self.ymin + num.arange(self.ny) * self.dy
85 def positive_region(self, region):
86 xmin, xmax, ymin, ymax = [float(x) for x in region]
88 assert -180. - 360. <= xmin < 180.
89 assert -180. < xmax <= 180. + 360.
90 assert -90. <= ymin < 90.
91 assert -90. < ymax <= 90.
93 if xmax < xmin:
94 xmax += 360.
96 if xmin < -180.:
97 xmin += 360.
98 xmax += 360.
100 return (xmin, xmax, ymin, ymax)
102 def tile_indices(self, region):
103 xmin, xmax, ymin, ymax = self.positive_region(region)
104 itxmin = int(math.floor((xmin - self.xmin) / self.stx))
105 itxmax = int(math.ceil((xmax - self.xmin) / self.stx))
106 if itxmin == itxmax:
107 itxmax += 1
108 itymin = int(math.floor((ymin - self.ymin) / self.sty))
109 itymax = int(math.ceil((ymax - self.ymin) / self.sty))
110 if itymin == itymax:
111 if itymax != self.ntilesy:
112 itymax += 1
113 else:
114 itymin -= 1
116 indices = []
117 for ity in range(itymin, itymax):
118 for itx in range(itxmin, itxmax):
119 indices.append((itx % self.ntilesx, ity))
121 return indices
123 def get_tile(self, itx, ity):
124 return None
126 def get(self, region):
127 if len(region) == 2:
128 x, y = region
129 t = self.get((x, x, y, y))
130 if t is not None:
131 return t.get(x, y)
132 else:
133 return None
135 indices = self.tile_indices(region)
136 tiles = []
137 for itx, ity in indices:
138 t = self.get_tile(itx, ity)
139 if t:
140 tiles.append(t)
142 return tile.combine(tiles, region)
144 def get_with_repeat(self, region):
145 xmin, xmax, ymin, ymax = region
146 ymin2 = max(-90., ymin)
147 ymax2 = min(90., ymax)
148 region2 = xmin, xmax, ymin2, ymax2
149 t = self.get(region2)
150 if t is not None and region2 != region:
151 t.yextend_with_repeat(ymin, ymax)
153 return t
156class DecimatedTiledGlobalDataset(TiledGlobalDataset):
158 def __init__(self, name, base, ndeci, data_dir=None, ntx=None, nty=None):
160 assert ndeci % 2 == 0
161 assert (base.nx - 1) % ndeci == 0
162 assert (base.ny - 1) % ndeci == 0
164 nx = (base.nx - 1) // ndeci + 1
165 ny = (base.ny - 1) // ndeci + 1
167 if ntx is None:
168 ntx = base.ntx
170 if nty is None:
171 nty = base.nty
173 assert (nx - 1) % (ntx - 1) == 0
174 assert (ny - 1) % (nty - 1) == 0
176 if data_dir is None:
177 data_dir = op.join(base.data_dir, 'decimated_%i' % ndeci)
179 TiledGlobalDataset.__init__(self, name, nx, ny, ntx, nty, base.dtype,
180 data_dir=data_dir, citation=base.citation,
181 region=base.region)
183 self.base = base
184 self.ndeci = ndeci
186 def make_tile(self, itx, ity, fpath):
187 nh = self.ndeci // 2
188 xmin = self.xmin + itx*self.stx - self.base.dx * nh
189 xmax = self.xmin + (itx+1)*self.stx + self.base.dx * nh
190 ymin = self.ymin + ity*self.sty - self.base.dy * nh
191 ymax = self.ymin + (ity+1)*self.sty + self.base.dy * nh
193 t = self.base.get_with_repeat((xmin, xmax, ymin, ymax))
194 if t is not None:
195 t.decimate(self.ndeci)
197 util.ensuredirs(fpath)
198 with open(fpath, 'w') as f:
199 if t is not None:
200 t.data.tofile(f)
202 def make_if_needed(self, itx, ity):
203 assert 0 <= itx < self.ntilesx
204 assert 0 <= ity < self.ntilesy
206 fn = '%02i.%02i.bin' % (ity, itx)
207 fpath = op.join(self.data_dir, fn)
208 if not op.exists(fpath):
209 logger.info('making decimated tile: %s (%s)' % (fn, self.name))
210 self.make_tile(itx, ity, fpath)
212 def get_tile(self, itx, ity):
213 assert 0 <= itx < self.ntilesx
214 assert 0 <= ity < self.ntilesy
216 self.make_if_needed(itx, ity)
218 fn = '%02i.%02i.bin' % (ity, itx)
219 fpath = op.join(self.data_dir, fn)
220 with open(fpath, 'r') as f:
221 data = num.fromfile(f, dtype=self.dtype)
223 if data.size == 0:
224 return None
226 assert data.size == self.ntx*self.nty
227 data = data.reshape((self.ntx, self.nty))
229 return tile.Tile(
230 self.xmin + itx*self.stx,
231 self.ymin + ity*self.sty,
232 self.dx, self.dx, data)
234 def make_all_missing(self):
235 for ity in range(self.ntilesy):
236 for itx in range(self.ntilesx):
237 self.make_if_needed(itx, ity)