1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import zipfile 

7import os.path as op 

8 

9import numpy as num 

10 

11from . import tile, dataset 

12 

13try: 

14 range = xrange 

15except NameError: 

16 pass 

17 

18citation = ''' 

19Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: 

20Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS 

21NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M 

22[access date]. 

23''' 

24 

25 

26class ETOPO1(dataset.TiledGlobalDataset): 

27 

28 def __init__( 

29 self, 

30 name='ETOPO1', 

31 data_dir=op.join(op.dirname(__file__), 'data', 'ETOPO1'), 

32 base_fn='etopo1_ice_g_i2', 

33 raw_data_url=('https://mirror.pyrocko.org/www.ngdc.noaa.gov/mgg' 

34 '/global/relief/ETOPO1/data/ice_surface/' 

35 'grid_registered/binary/%s.zip')): 

36 

37 dataset.TiledGlobalDataset.__init__( 

38 self, 

39 name, 

40 21601, 10801, 1351, 1351, 

41 num.dtype('<i2'), 

42 data_dir=data_dir, 

43 citation=citation) 

44 

45 self.base_fn = base_fn 

46 self.raw_data_url = raw_data_url 

47 

48 def download(self): 

49 fpath = op.join(self.data_dir, '%s.zip' % self.base_fn) 

50 if not op.exists(fpath): 

51 self.download_file( 

52 self.raw_data_url % self.base_fn, fpath) 

53 

54 self.make_tiles() 

55 

56 def make_tiles(self): 

57 fpath = op.join(self.data_dir, '%s.zip' % self.base_fn) 

58 

59 zipf = zipfile.ZipFile(fpath, 'r') 

60 rawdata = zipf.read('%s.bin' % self.base_fn) 

61 zipf.close() 

62 data = num.frombuffer(rawdata, dtype=self.dtype) 

63 assert data.size == self.nx * self.ny 

64 data = data.reshape((self.ny, self.nx))[::-1] 

65 

66 for ity in range(self.ntilesy): 

67 for itx in range(self.ntilesx): 

68 tiledata = data[ity*(self.nty-1):(ity+1)*(self.nty-1)+1, 

69 itx*(self.ntx-1):(itx+1)*(self.ntx-1)+1] 

70 

71 fn = '%s.%02i.%02i.bin' % (self.base_fn, ity, itx) 

72 fpath = op.join(self.data_dir, fn) 

73 with open(fpath, 'w') as f: 

74 tiledata.tofile(f) 

75 

76 def get_tile(self, itx, ity): 

77 assert 0 <= itx < self.ntilesx 

78 assert 0 <= ity < self.ntilesy 

79 

80 fn = '%s.%02i.%02i.bin' % (self.base_fn, ity, itx) 

81 fpath = op.join(self.data_dir, fn) 

82 if not op.exists(fpath): 

83 self.download() 

84 

85 with open(fpath, 'r') as f: 

86 data = num.fromfile(f, dtype=self.dtype) 

87 

88 assert data.size == self.ntx*self.nty 

89 data = data.reshape((self.ntx, self.nty)) 

90 

91 return tile.Tile( 

92 self.xmin + itx*self.stx, 

93 self.ymin + ity*self.sty, 

94 self.dx, self.dx, data)