Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/dataset/topo/etopo1.py: 40%

48 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7NOAA `ETOPO1 <https://www.ngdc.noaa.gov/mgg/global/>`_ Global Relief Map. 

8''' 

9 

10import zipfile 

11import os.path as op 

12 

13import numpy as num 

14 

15from . import tile, dataset 

16 

17try: 

18 range = xrange 

19except NameError: 

20 pass 

21 

22citation = ''' 

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

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

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

26''' 

27 

28__doc__ += ''' 

29 

30.. rubric:: Citation 

31''' + citation 

32 

33 

34class ETOPO1(dataset.TiledGlobalDataset): 

35 ''' 

36 ETOPO1 as a tiled global dataset. 

37 ''' 

38 

39 def __init__( 

40 self, 

41 name='ETOPO1', 

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

43 base_fn='etopo1_ice_g_i2', 

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

45 '/global/relief/ETOPO1/data/ice_surface/' 

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

47 

48 dataset.TiledGlobalDataset.__init__( 

49 self, 

50 name, 

51 21601, 10801, 1351, 1351, 

52 num.dtype('<i2'), 

53 data_dir=data_dir, 

54 citation=citation) 

55 

56 self.base_fn = base_fn 

57 self.raw_data_url = raw_data_url 

58 

59 def download(self): 

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

61 if not op.exists(fpath): 

62 self.download_file( 

63 self.raw_data_url % self.base_fn, fpath) 

64 

65 self.make_tiles() 

66 

67 def make_tiles(self): 

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

69 

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

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

72 zipf.close() 

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

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

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

76 

77 for ity in range(self.ntilesy): 

78 for itx in range(self.ntilesx): 

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

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

81 

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

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

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

85 tiledata.tofile(f) 

86 

87 def get_tile(self, itx, ity): 

88 assert 0 <= itx < self.ntilesx 

89 assert 0 <= ity < self.ntilesy 

90 

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

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

93 if not op.exists(fpath): 

94 self.download() 

95 

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

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

98 

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

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

101 

102 return tile.Tile( 

103 self.xmin + itx*self.stx, 

104 self.ymin + ity*self.sty, 

105 self.dx, self.dx, data)