1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import
7import zipfile
8import os.path as op
9import os
10import re
11# from pyrocko.util import urlopen
13import numpy as num
15from pyrocko import util, config
16from . import tile, dataset
18citation = '''
19Farr, T. G., and M. Kobrick, 2000, Shuttle Radar Topography Mission produces a
20wealth of data. Eos Trans. AGU, 81:583-585.
22Farr, T. G. et al., 2007, The Shuttle Radar Topography Mission, Rev. Geophys.,
2345, RG2004, doi:10.1029/2005RG000183. (Also available online at
24http://www2.jpl.nasa.gov/srtm/SRTM_paper.pdf)
26Kobrick, M., 2006, On the toes of giants-How SRTM was born, Photogramm. Eng.
27Remote Sens., 72:206-210.
29Rosen, P. A. et al., 2000, Synthetic aperture radar interferometry, Proc. IEEE,
3088:333-382.
31'''
34class AuthenticationRequired(Exception):
35 pass
38class SRTMGL3(dataset.TiledGlobalDataset):
40 def __init__(
41 self,
42 name='SRTMGL3',
43 data_dir=op.join(op.dirname(__file__), 'data', 'SRTMGL3'),
44 raw_data_url='https://mirror.pyrocko.org/e4ftl01.cr.usgs.gov/'
45 'MEASURES/SRTMGL3.003/2000.02.11'):
47 dataset.TiledGlobalDataset.__init__(
48 self,
49 name,
50 432001, 216001, 1201, 1201,
51 num.dtype('>i2'),
52 data_dir=data_dir,
53 citation=citation,
54 region=(-180., 180., -60., 60.))
56 self.raw_data_url = raw_data_url
57 self._available_tilenames = None
58 self.config = config.config()
60 def tilename(self, itx, ity):
61 itx -= 180
62 ity -= 90
63 if ity >= 0:
64 s = 'N%02i' % ity
65 else:
66 s = 'S%02i' % -ity
68 if itx >= 0:
69 s += 'E%03i' % itx
70 else:
71 s += 'W%03i' % -itx
73 return s
75 def available_tilenames(self):
76 if self._available_tilenames is None:
77 fpath = op.join(self.data_dir, 'available.list')
78 if not op.exists(fpath) or os.stat(fpath).st_size == 0:
79 util.ensuredirs(fpath)
80 # remote structure changed, we would have to clawl through
81 # many pages. Now keeping tile index here:
82 self.download_file(
83 'https://mirror.pyrocko.org/e4ftl01.cr.usgs.gov/'
84 'MEASURES/SRTMGL3.003/2000.02.11/available.list', fpath)
86 # url = self.raw_data_url + '/'
87 # f = urlopen(url)
88 # data = f.read().decode()
89 # available = re.findall(
90 # r'([NS]\d\d[EW]\d\d\d)\.SRTMGL3\.hgt', data)
91 #
92 # f.close()
93 #
94 # with open(fpath, 'w') as f:
95 # f.writelines('%s\n' % s for s in available)
97 with open(fpath, 'r') as f:
98 available = [
99 s.strip() for s in f.readlines()
100 if re.match(r'^[NS]\d\d[EW]\d\d\d$', s.strip())]
102 self._available_tilenames = set(available)
104 return self._available_tilenames
106 def tilefilename(self, tilename):
107 return tilename + '.SRTMGL3.hgt.zip'
109 def tilepath(self, tilename):
110 fn = self.tilefilename(tilename)
111 return op.join(self.data_dir, fn)
113 def download_tile(self, tilename):
114 fpath = self.tilepath(tilename)
115 fn = self.tilefilename(tilename)
116 url = self.raw_data_url + '/' + fn
117 try:
118 # we have to follow the oauth redirect here...
119 self.download_file(url, fpath)
120 except Exception as e:
121 raise e
123 def download(self):
124 for tn in self.available_tilenames():
125 fpath = self.tilepath(tn)
126 if not op.exists(fpath):
127 self.download_tile(tn)
129 def get_tile(self, itx, ity):
130 tn = self.tilename(itx, ity)
131 if tn not in self.available_tilenames():
132 return None
133 else:
134 fpath = self.tilepath(tn)
135 if not op.exists(fpath):
136 self.download_tile(tn)
138 zipf = zipfile.ZipFile(fpath, 'r')
139 rawdata = zipf.read(tn + '.hgt')
140 zipf.close()
141 data = num.fromstring(rawdata, dtype=self.dtype)
142 assert data.size == self.ntx * self.nty
143 data = data.reshape(self.nty, self.ntx)[::-1, ::]
144 return tile.Tile(
145 self.xmin + itx*self.stx,
146 self.ymin + ity*self.sty,
147 self.dx, self.dx, data)
150if __name__ == '__main__':
151 import sys
153 util.setup_logging('pyrocko.topo.srtmgl3', 'info')
155 if len(sys.argv) != 2:
156 sys.exit('usage: python -m pyrocko.topo.srtmgl3 download')
158 if sys.argv[1] == 'download':
159 srtmgl3 = SRTMGL3()
160 srtmgl3.download()