• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

OGGM / oggm / 3836586642

pending completion
3836586642

Pull #1510

github

GitHub
Merge 2a19b5c72 into 8ab12c504
Pull Request #1510: Fix issue with TANDEM DEM and minor RGI7 issue

2 of 4 new or added lines in 2 files covered. (50.0%)

249 existing lines in 7 files now uncovered.

12099 of 13681 relevant lines covered (88.44%)

3.61 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

78.09
/oggm/utils/_downloads.py
1
"""Automated data download and IO."""
2

3
# Builtins
4
import glob
9✔
5
import os
9✔
6
import gzip
9✔
7
import bz2
9✔
8
import hashlib
9✔
9
import shutil
9✔
10
import zipfile
9✔
11
import sys
9✔
12
import math
9✔
13
import logging
9✔
14
from functools import partial, wraps
9✔
15
import time
9✔
16
import fnmatch
9✔
17
import urllib.request
9✔
18
import urllib.error
9✔
19
from urllib.parse import urlparse
9✔
20
import socket
9✔
21
import multiprocessing
9✔
22
from netrc import netrc
9✔
23
import ftplib
9✔
24
import ssl
9✔
25
import tarfile
9✔
26

27
# External libs
28
import pandas as pd
9✔
29
import numpy as np
9✔
30
import shapely.geometry as shpg
9✔
31
import requests
9✔
32

33
# Optional libs
34
try:
9✔
35
    import geopandas as gpd
9✔
36
except ImportError:
37
    pass
38
try:
9✔
39
    import salem
9✔
40
    from salem import wgs84
8✔
41
except ImportError:
42
    pass
43
try:
9✔
44
    import rasterio
9✔
45
    try:
8✔
46
        # rasterio V > 1.0
47
        from rasterio.merge import merge as merge_tool
8✔
48
    except ImportError:
49
        from rasterio.tools.merge import merge as merge_tool
50
except ImportError:
51
    pass
52
try:
9✔
53
    ModuleNotFoundError
9✔
54
except NameError:
×
55
    ModuleNotFoundError = ImportError
×
56

57
# Locals
58
import oggm.cfg as cfg
9✔
59
from oggm.exceptions import (InvalidParamsError, NoInternetException,
9✔
60
                             DownloadVerificationFailedException,
61
                             DownloadCredentialsMissingException,
62
                             HttpDownloadError, HttpContentTooShortError,
63
                             InvalidDEMError, FTPSDownloadError)
64

65
# Module logger
66
logger = logging.getLogger('.'.join(__name__.split('.')[:-1]))
9✔
67

68
# Github repository and commit hash/branch name/tag name on that repository
69
# The given commit will be downloaded from github and used as source for
70
# all sample data
71
SAMPLE_DATA_GH_REPO = 'OGGM/oggm-sample-data'
9✔
72
SAMPLE_DATA_COMMIT = 'd9f01846960a141690bab9d38a85524d89a0d9ae'
9✔
73

74
CHECKSUM_URL = 'https://cluster.klima.uni-bremen.de/data/downloads.sha256.hdf'
9✔
75
CHECKSUM_VALIDATION_URL = CHECKSUM_URL + '.sha256'
9✔
76
CHECKSUM_LIFETIME = 24 * 60 * 60
9✔
77

78
# Web mercator proj constants
79
WEB_N_PIX = 256
9✔
80
WEB_EARTH_RADUIS = 6378137.
9✔
81

82
DEM_SOURCES = ['GIMP', 'ARCTICDEM', 'RAMP', 'TANDEM', 'AW3D30', 'MAPZEN', 'DEM3',
9✔
83
               'ASTER', 'SRTM', 'REMA', 'ALASKA', 'COPDEM30', 'COPDEM90', 'NASADEM']
84
DEM_SOURCES_PER_GLACIER = None
9✔
85

86
_RGI_METADATA = dict()
9✔
87

88
DEM3REG = {
9✔
89
    'ISL': [-25., -13., 63., 67.],  # Iceland
90
    'SVALBARD': [9., 35.99, 75., 84.],
91
    'JANMAYEN': [-10., -7., 70., 72.],
92
    'FJ': [36., 68., 79., 90.],  # Franz Josef Land
93
    'FAR': [-8., -6., 61., 63.],  # Faroer
94
    'BEAR': [18., 20., 74., 75.],  # Bear Island
95
    'SHL': [-3., 0., 60., 61.],  # Shetland
96
    # Antarctica tiles as UTM zones, large files
97
    '01-15': [-180., -91., -90, -60.],
98
    '16-30': [-91., -1., -90., -60.],
99
    '31-45': [-1., 89., -90., -60.],
100
    '46-60': [89., 189., -90., -60.],
101
    # Greenland tiles
102
    'GL-North': [-72., -11., 76., 84.],
103
    'GL-West': [-62., -42., 64., 76.],
104
    'GL-South': [-52., -40., 59., 64.],
105
    'GL-East': [-42., -17., 64., 76.]
106
}
107

108
# Function
109
tuple2int = partial(np.array, dtype=np.int64)
9✔
110

111
lock = None
9✔
112

113

114
def mkdir(path, reset=False):
9✔
115
    """Checks if directory exists and if not, create one.
116

117
    Parameters
118
    ----------
119
    reset: erase the content of the directory if exists
120

121
    Returns
122
    -------
123
    the path
124
    """
125

126
    if reset and os.path.exists(path):
8✔
127
        shutil.rmtree(path)
×
128
    try:
8✔
129
        os.makedirs(path)
8✔
130
    except FileExistsError:
8✔
131
        pass
8✔
132
    return path
8✔
133

134

135
def del_empty_dirs(s_dir):
9✔
136
    """Delete empty directories."""
137
    b_empty = True
1✔
138
    for s_target in os.listdir(s_dir):
1✔
139
        s_path = os.path.join(s_dir, s_target)
1✔
140
        if os.path.isdir(s_path):
1✔
141
            if not del_empty_dirs(s_path):
1✔
142
                b_empty = False
×
143
        else:
144
            b_empty = False
1✔
145
    if b_empty:
1✔
146
        os.rmdir(s_dir)
1✔
147
    return b_empty
1✔
148

149

150
def findfiles(root_dir, endswith):
9✔
151
    """Finds all files with a specific ending in a directory
152

153
    Parameters
154
    ----------
155
    root_dir : str
156
       The directory to search fo
157
    endswith : str
158
       The file ending (e.g. '.hgt'
159

160
    Returns
161
    -------
162
    the list of files
163
    """
164
    out = []
×
165
    for dirpath, dirnames, filenames in os.walk(root_dir):
×
166
        for filename in [f for f in filenames if f.endswith(endswith)]:
×
167
            out.append(os.path.join(dirpath, filename))
×
168
    return out
×
169

170

171
def get_lock():
9✔
172
    """Get multiprocessing lock."""
173
    global lock
174
    if lock is None:
8✔
175
        # Global Lock
176
        if cfg.PARAMS.get('use_mp_spawn', False):
8✔
177
            lock = multiprocessing.get_context('spawn').Lock()
×
178
        else:
179
            lock = multiprocessing.Lock()
8✔
180
    return lock
8✔
181

182

183
def get_dl_verify_data(section):
9✔
184
    """Returns a pandas DataFrame with all known download object hashes.
185

186
    The returned dictionary resolves str: cache_obj_name (without section)
187
    to a tuple int(size) and bytes(sha256)
188
    """
189

190
    verify_key = 'dl_verify_data_' + section
8✔
191
    if verify_key in cfg.DATA:
8✔
192
        return cfg.DATA[verify_key]
4✔
193

194
    verify_file_path = os.path.join(cfg.CACHE_DIR, 'downloads.sha256.hdf')
8✔
195

196
    def verify_file(force=False):
8✔
197
        """Check the hash file's own hash"""
198
        if not cfg.PARAMS['has_internet']:
8✔
199
            return
×
200

201
        if not force and os.path.isfile(verify_file_path) and \
8✔
202
           os.path.getmtime(verify_file_path) + CHECKSUM_LIFETIME > time.time():
203
            return
7✔
204

205
        logger.info('Checking the download verification file checksum...')
8✔
206
        try:
8✔
207
            with requests.get(CHECKSUM_VALIDATION_URL) as req:
8✔
208
                req.raise_for_status()
8✔
209
                verify_file_sha256 = req.text.split(maxsplit=1)[0]
8✔
210
                verify_file_sha256 = bytearray.fromhex(verify_file_sha256)
8✔
211
        except Exception as e:
×
212
            verify_file_sha256 = None
×
213
            logger.warning('Failed getting verification checksum: ' + repr(e))
×
214

215
        if os.path.isfile(verify_file_path) and verify_file_sha256:
8✔
216
            sha256 = hashlib.sha256()
8✔
217
            with open(verify_file_path, 'rb') as f:
8✔
218
                for b in iter(lambda: f.read(0xFFFF), b''):
8✔
219
                    sha256.update(b)
8✔
220
            if sha256.digest() != verify_file_sha256:
8✔
221
                logger.warning('%s changed or invalid, deleting.'
×
222
                               % (verify_file_path))
223
                os.remove(verify_file_path)
×
224
            else:
225
                os.utime(verify_file_path)
8✔
226

227
    if not np.any(['dl_verify_data_' in k for k in cfg.DATA.keys()]):
8✔
228
        # We check the hash file only once per session
229
        # no need to do it at each call
230
        verify_file()
8✔
231

232
    if not os.path.isfile(verify_file_path):
8✔
233
        if not cfg.PARAMS['has_internet']:
8✔
234
            return pd.DataFrame()
×
235

236
        logger.info('Downloading %s to %s...'
8✔
237
                    % (CHECKSUM_URL, verify_file_path))
238

239
        with requests.get(CHECKSUM_URL, stream=True) as req:
8✔
240
            if req.status_code == 200:
8✔
241
                mkdir(os.path.dirname(verify_file_path))
8✔
242
                with open(verify_file_path, 'wb') as f:
8✔
243
                    for b in req.iter_content(chunk_size=0xFFFF):
8✔
244
                        if b:
8✔
245
                            f.write(b)
8✔
246

247
        logger.info('Done downloading.')
8✔
248

249
        verify_file(force=True)
8✔
250

251
    if not os.path.isfile(verify_file_path):
8✔
252
        logger.warning('Downloading and verifying checksums failed.')
×
253
        return pd.DataFrame()
×
254

255
    try:
8✔
256
        data = pd.read_hdf(verify_file_path, key=section)
8✔
257
    except KeyError:
8✔
258
        data = pd.DataFrame()
8✔
259

260
    cfg.DATA[verify_key] = data
8✔
261

262
    return data
8✔
263

264

265
def _call_dl_func(dl_func, cache_path):
9✔
266
    """Helper so the actual call to downloads can be overridden
267
    """
268
    return dl_func(cache_path)
8✔
269

270

271
def _cached_download_helper(cache_obj_name, dl_func, reset=False):
9✔
272
    """Helper function for downloads.
273

274
    Takes care of checking if the file is already cached.
275
    Only calls the actual download function when no cached version exists.
276
    """
277
    cache_dir = cfg.PATHS['dl_cache_dir']
8✔
278
    cache_ro = cfg.PARAMS['dl_cache_readonly']
8✔
279

280
    # A lot of logic below could be simplified but it's also not too important
281
    wd = cfg.PATHS.get('working_dir')
8✔
282
    if wd:
8✔
283
        # this is for real runs
284
        fb_cache_dir = os.path.join(wd, 'cache')
7✔
285
        check_fb_dir = False
7✔
286
    else:
287
        # Nothing have been set up yet, this is bad - find a place to write
288
        # This should happen on read-only cluster only but still
289
        wd = os.environ.get('OGGM_WORKDIR')
8✔
290
        if wd is not None and os.path.isdir(wd):
8✔
291
            fb_cache_dir = os.path.join(wd, 'cache')
×
292
        else:
293
            fb_cache_dir = os.path.join(cfg.CACHE_DIR, 'cache')
8✔
294
        check_fb_dir = True
8✔
295

296
    if not cache_dir:
8✔
297
        # Defaults to working directory: it must be set!
298
        if not cfg.PATHS['working_dir']:
×
299
            raise InvalidParamsError("Need a valid PATHS['working_dir']!")
×
300
        cache_dir = fb_cache_dir
×
301
        cache_ro = False
×
302

303
    fb_path = os.path.join(fb_cache_dir, cache_obj_name)
8✔
304
    if not reset and os.path.isfile(fb_path):
8✔
305
        return fb_path
×
306

307
    cache_path = os.path.join(cache_dir, cache_obj_name)
8✔
308
    if not reset and os.path.isfile(cache_path):
8✔
309
        return cache_path
6✔
310

311
    if cache_ro:
8✔
312
        if check_fb_dir:
×
313
            # Add a manual check that we are caching sample data download
314
            if 'oggm-sample-data' not in fb_path:
×
315
                raise InvalidParamsError('Attempting to download something '
×
316
                                         'with invalid global settings.')
317
        cache_path = fb_path
×
318

319
    if not cfg.PARAMS['has_internet']:
8✔
320
        raise NoInternetException("Download required, but "
1✔
321
                                  "`has_internet` is False.")
322

323
    mkdir(os.path.dirname(cache_path))
8✔
324

325
    try:
8✔
326
        cache_path = _call_dl_func(dl_func, cache_path)
8✔
327
    except BaseException:
×
328
        if os.path.exists(cache_path):
×
329
            os.remove(cache_path)
×
330
        raise
×
331

332
    return cache_path
8✔
333

334

335
def _verified_download_helper(cache_obj_name, dl_func, reset=False):
9✔
336
    """Helper function for downloads.
337

338
    Verifies the size and hash of the downloaded file against the included
339
    list of known static files.
340
    Uses _cached_download_helper to perform the actual download.
341
    """
342
    path = _cached_download_helper(cache_obj_name, dl_func, reset)
8✔
343

344
    try:
8✔
345
        dl_verify = cfg.PARAMS['dl_verify']
8✔
346
    except KeyError:
1✔
347
        dl_verify = True
1✔
348

349
    if dl_verify and path and cache_obj_name not in cfg.DL_VERIFIED:
8✔
350
        cache_section, cache_path = cache_obj_name.split('/', 1)
8✔
351
        data = get_dl_verify_data(cache_section)
8✔
352
        if cache_path not in data.index:
8✔
353
            logger.info('No known hash for %s' % cache_obj_name)
8✔
354
            cfg.DL_VERIFIED[cache_obj_name] = True
8✔
355
        else:
356
            # compute the hash
357
            sha256 = hashlib.sha256()
7✔
358
            with open(path, 'rb') as f:
7✔
359
                for b in iter(lambda: f.read(0xFFFF), b''):
7✔
360
                    sha256.update(b)
7✔
361
            sha256 = sha256.digest()
7✔
362
            size = os.path.getsize(path)
7✔
363

364
            # check
365
            data = data.loc[cache_path]
7✔
366
            if data['size'] != size or bytes(data['sha256']) != sha256:
7✔
367
                err = '%s failed to verify!\nis: %s %s\nexpected: %s %s' % (
1✔
368
                    path, size, sha256.hex(), data[0], data[1].hex())
369
                raise DownloadVerificationFailedException(msg=err, path=path)
1✔
370
            logger.info('%s verified successfully.' % path)
7✔
371
            cfg.DL_VERIFIED[cache_obj_name] = True
7✔
372

373
    return path
8✔
374

375

376
def _requests_urlretrieve(url, path, reporthook, auth=None, timeout=None):
9✔
377
    """Implements the required features of urlretrieve on top of requests
378
    """
379

380
    chunk_size = 128 * 1024
8✔
381
    chunk_count = 0
8✔
382

383
    with requests.get(url, stream=True, auth=auth, timeout=timeout) as r:
8✔
384
        if r.status_code != 200:
8✔
385
            raise HttpDownloadError(r.status_code, url)
×
386
        r.raise_for_status()
8✔
387

388
        size = r.headers.get('content-length') or -1
8✔
389
        size = int(size)
8✔
390

391
        if reporthook:
8✔
392
            reporthook(chunk_count, chunk_size, size)
8✔
393

394
        with open(path, 'wb') as f:
8✔
395
            for chunk in r.iter_content(chunk_size=chunk_size):
8✔
396
                if not chunk:
8✔
397
                    continue
×
398
                f.write(chunk)
8✔
399
                chunk_count += 1
8✔
400
                if reporthook:
8✔
401
                    reporthook(chunk_count, chunk_size, size)
8✔
402

403
        if chunk_count * chunk_size < size:
8✔
404
            raise HttpContentTooShortError()
×
405

406

407
def _classic_urlretrieve(url, path, reporthook, auth=None, timeout=None):
9✔
408
    """Thin wrapper around pythons urllib urlretrieve
409
    """
410

411
    ourl = url
×
412
    if auth:
×
413
        u = urlparse(url)
×
414
        if '@' not in u.netloc:
×
415
            netloc = auth[0] + ':' + auth[1] + '@' + u.netloc
×
416
            url = u._replace(netloc=netloc).geturl()
×
417

418
    old_def_timeout = socket.getdefaulttimeout()
×
419
    if timeout is not None:
×
420
        socket.setdefaulttimeout(timeout)
×
421

422
    try:
×
423
        urllib.request.urlretrieve(url, path, reporthook)
×
424
    except urllib.error.HTTPError as e:
×
425
        raise HttpDownloadError(e.code, ourl)
×
426
    except urllib.error.ContentTooShortError as e:
×
427
        raise HttpContentTooShortError()
×
428
    finally:
429
        socket.setdefaulttimeout(old_def_timeout)
×
430

431

432
class ImplicitFTPTLS(ftplib.FTP_TLS):
9✔
433
    """ FTP_TLS subclass that automatically wraps sockets in SSL to support
434
        implicit FTPS.
435

436
        Taken from https://stackoverflow.com/a/36049814
437
    """
438

439
    def __init__(self, *args, **kwargs):
9✔
440
        super().__init__(*args, **kwargs)
×
441
        self._sock = None
×
442

443
    @property
9✔
444
    def sock(self):
9✔
445
        """Return the socket."""
446
        return self._sock
×
447

448
    @sock.setter
9✔
449
    def sock(self, value):
9✔
450
        """When modifying the socket, ensure that it is ssl wrapped."""
451
        if value is not None and not isinstance(value, ssl.SSLSocket):
×
452
            value = self.context.wrap_socket(value)
×
453
        self._sock = value
×
454

455

456
def url_exists(url):
9✔
457
    """Checks if a given a URL exists or not."""
458
    request = requests.get(url)
2✔
459
    return request.status_code < 400
2✔
460

461

462
def _ftps_retrieve(url, path, reporthook, auth=None, timeout=None):
9✔
463
    """ Wrapper around ftplib to download from FTPS server
464
    """
465

466
    if not auth:
×
467
        raise DownloadCredentialsMissingException('No authentication '
×
468
                                                  'credentials given!')
469

470
    upar = urlparse(url)
×
471

472
    # Decide if Implicit or Explicit FTPS is used based on the port in url
473
    if upar.port == 990:
×
474
        ftps = ImplicitFTPTLS()
×
475
    elif upar.port == 21:
×
476
        ftps = ftplib.FTP_TLS()
×
477

478
    try:
×
479
        # establish ssl connection
480
        ftps.connect(host=upar.hostname, port=upar.port, timeout=timeout)
×
481
        ftps.login(user=auth[0], passwd=auth[1])
×
482
        ftps.prot_p()
×
483

484
        logger.info('Established connection %s' % upar.hostname)
×
485

486
        # meta for progress bar size
487
        count = 0
×
488
        total = ftps.size(upar.path)
×
489
        bs = 12*1024
×
490

491
        def _ftps_progress(data):
×
492
            outfile.write(data)
×
493
            nonlocal count
494
            count += 1
×
495
            reporthook(count, count*bs, total)
×
496

497
        with open(path, 'wb') as outfile:
×
498
            ftps.retrbinary('RETR ' + upar.path, _ftps_progress, blocksize=bs)
×
499

500
    except (ftplib.error_perm, socket.timeout, socket.gaierror) as err:
×
501
        raise FTPSDownloadError(err)
×
502
    finally:
503
        ftps.close()
×
504

505

506
def _get_url_cache_name(url):
9✔
507
    """Returns the cache name for any given url.
508
    """
509

510
    res = urlparse(url)
8✔
511
    return res.netloc.split(':', 1)[0] + res.path
8✔
512

513

514
def oggm_urlretrieve(url, cache_obj_name=None, reset=False,
9✔
515
                     reporthook=None, auth=None, timeout=None):
516
    """Wrapper around urlretrieve, to implement our caching logic.
517

518
    Instead of accepting a destination path, it decided where to store the file
519
    and returns the local path.
520

521
    auth is expected to be either a tuple of ('username', 'password') or None.
522
    """
523

524
    if cache_obj_name is None:
8✔
525
        cache_obj_name = _get_url_cache_name(url)
8✔
526

527
    def _dlf(cache_path):
8✔
528
        logger.info("Downloading %s to %s..." % (url, cache_path))
8✔
529
        try:
8✔
530
            _requests_urlretrieve(url, cache_path, reporthook, auth, timeout)
8✔
531
        except requests.exceptions.InvalidSchema:
×
532
            if 'ftps://' in url:
×
533
                _ftps_retrieve(url, cache_path, reporthook, auth, timeout)
×
534
            else:
535
                _classic_urlretrieve(url, cache_path, reporthook, auth,
×
536
                                     timeout)
537
        return cache_path
8✔
538

539
    return _verified_download_helper(cache_obj_name, _dlf, reset)
8✔
540

541

542
def _progress_urlretrieve(url, cache_name=None, reset=False,
9✔
543
                          auth=None, timeout=None):
544
    """Downloads a file, returns its local path, and shows a progressbar."""
545

546
    try:
8✔
547
        from progressbar import DataTransferBar, UnknownLength
8✔
548
        pbar = None
8✔
549

550
        def _upd(count, size, total):
8✔
551
            nonlocal pbar
552
            if pbar is None:
8✔
553
                pbar = DataTransferBar()
8✔
554
                if not pbar.is_terminal:
8✔
555
                    pbar.min_poll_interval = 15
8✔
556
            if pbar.max_value is None:
8✔
557
                if total > 0:
8✔
558
                    pbar.start(total)
8✔
559
                else:
UNCOV
560
                    pbar.start(UnknownLength)
×
561
            pbar.update(min(count * size, total))
8✔
562
            sys.stdout.flush()
8✔
563
        res = oggm_urlretrieve(url, cache_obj_name=cache_name, reset=reset,
8✔
564
                               reporthook=_upd, auth=auth, timeout=timeout)
565
        try:
8✔
566
            pbar.finish()
8✔
567
        except BaseException:
6✔
568
            pass
6✔
569
        return res
8✔
570
    except (ImportError, ModuleNotFoundError):
1✔
571
        return oggm_urlretrieve(url, cache_obj_name=cache_name,
×
572
                                reset=reset, auth=auth, timeout=timeout)
573

574

575
def aws_file_download(aws_path, cache_name=None, reset=False):
9✔
576
    with get_lock():
×
577
        return _aws_file_download_unlocked(aws_path, cache_name, reset)
×
578

579

580
def _aws_file_download_unlocked(aws_path, cache_name=None, reset=False):
9✔
581
    """Download a file from the AWS drive s3://astgtmv2/
582

583
    **Note:** you need AWS credentials for this to work.
584

585
    Parameters
586
    ----------
587
    aws_path: path relative to s3://astgtmv2/
588
    """
589

590
    while aws_path.startswith('/'):
×
591
        aws_path = aws_path[1:]
×
592

593
    if cache_name is not None:
×
594
        cache_obj_name = cache_name
×
595
    else:
596
        cache_obj_name = 'astgtmv2/' + aws_path
×
597

598
    def _dlf(cache_path):
×
599
        raise NotImplementedError("Downloads from AWS are no longer supported")
600

601
    return _verified_download_helper(cache_obj_name, _dlf, reset)
×
602

603

604
def file_downloader(www_path, retry_max=3, sleep_on_retry=5,
9✔
605
                    cache_name=None, reset=False, auth=None,
606
                    timeout=None):
607
    """A slightly better downloader: it tries more than once."""
608

609
    local_path = None
8✔
610
    retry_counter = 0
8✔
611
    while retry_counter < retry_max:
8✔
612
        # Try to download
613
        try:
8✔
614
            retry_counter += 1
8✔
615
            local_path = _progress_urlretrieve(www_path, cache_name=cache_name,
8✔
616
                                               reset=reset, auth=auth,
617
                                               timeout=timeout)
618
            # if no error, exit
619
            break
8✔
620
        except HttpDownloadError as err:
1✔
621
            # This works well for py3
622
            if err.code == 404 or err.code == 300:
×
623
                # Ok so this *should* be an ocean tile
624
                return None
×
625
            elif err.code >= 500 and err.code < 600:
×
626
                logger.info(f"Downloading {www_path} failed with "
×
627
                            f"HTTP error {err.code}, "
628
                            f"retrying in {sleep_on_retry} seconds... "
629
                            f"{retry_counter}/{retry_max}")
630
                time.sleep(sleep_on_retry)
×
631
                continue
×
632
            else:
633
                raise
×
634
        except HttpContentTooShortError as err:
1✔
635
            logger.info("Downloading %s failed with ContentTooShortError"
×
636
                        " error %s, retrying in %s seconds... %s/%s" %
637
                        (www_path, err.code, sleep_on_retry, retry_counter, retry_max))
638
            time.sleep(sleep_on_retry)
×
639
            continue
×
640
        except DownloadVerificationFailedException as err:
1✔
641
            if (cfg.PATHS['dl_cache_dir'] and
1✔
642
                  err.path.startswith(cfg.PATHS['dl_cache_dir']) and
643
                  cfg.PARAMS['dl_cache_readonly']):
644
                if not cache_name:
×
645
                    cache_name = _get_url_cache_name(www_path)
×
646
                cache_name = "GLOBAL_CACHE_INVALID/" + cache_name
×
647
                retry_counter -= 1
×
648
                logger.info("Global cache for %s is invalid!")
×
649
            else:
650
                try:
1✔
651
                    os.remove(err.path)
1✔
652
                except FileNotFoundError:
×
653
                    pass
×
654
                logger.info("Downloading %s failed with "
1✔
655
                            "DownloadVerificationFailedException\n %s\n"
656
                            "The file might have changed or is corrupted. "
657
                            "File deleted. Re-downloading... %s/%s" %
658
                            (www_path, err.msg, retry_counter, retry_max))
659
            continue
1✔
660
        except requests.ConnectionError as err:
1✔
661
            if err.args[0].__class__.__name__ == 'MaxRetryError':
×
662
                # if request tried often enough we don't have to do this
663
                # this error does happen for not existing ASTERv3 files
664
                return None
×
665
            else:
666
                # in other cases: try again
667
                logger.info("Downloading %s failed with ConnectionError, "
×
668
                            "retrying in %s seconds... %s/%s" %
669
                            (www_path, sleep_on_retry, retry_counter, retry_max))
670
                time.sleep(sleep_on_retry)
×
671
                continue
×
672
        except FTPSDownloadError as err:
1✔
673
            logger.info("Downloading %s failed with FTPSDownloadError"
×
674
                        " error: '%s', retrying in %s seconds... %s/%s" %
675
                        (www_path, err.orgerr, sleep_on_retry, retry_counter, retry_max))
676
            time.sleep(sleep_on_retry)
×
677
            continue
×
678

679
    # See if we managed (fail is allowed)
680
    if not local_path or not os.path.exists(local_path):
8✔
681
        logger.warning('Downloading %s failed.' % www_path)
×
682

683
    return local_path
8✔
684

685

686
def locked_func(func):
9✔
687
    """To decorate a function that needs to be locked for multiprocessing"""
688
    @wraps(func)
9✔
689
    def wrapper(*args, **kwargs):
9✔
690
        with get_lock():
7✔
691
            return func(*args, **kwargs)
7✔
692
    return wrapper
9✔
693

694

695
def file_extractor(file_path):
9✔
696
    """For archives with only one file inside extract the file to tmpdir."""
697

698
    filename, file_extension = os.path.splitext(file_path)
7✔
699
    # Second one for tar.gz files
700
    f2, ex2 = os.path.splitext(filename)
7✔
701
    if ex2 == '.tar':
7✔
702
        filename, file_extension = f2, '.tar.gz'
1✔
703
    bname = os.path.basename(file_path)
7✔
704

705
    # This is to give a unique name to the tmp file
706
    hid = hashlib.md5(file_path.encode()).hexdigest()[:7] + '_'
7✔
707

708
    # extract directory
709
    tmpdir = cfg.PATHS['tmp_dir']
7✔
710
    mkdir(tmpdir)
7✔
711

712
    # Check output extension
713
    def _check_ext(f):
7✔
714
        _, of_ext = os.path.splitext(f)
7✔
715
        if of_ext not in ['.nc', '.tif']:
7✔
716
            raise InvalidParamsError('Extracted file extension not recognized'
×
717
                                     ': {}'.format(of_ext))
718
        return of_ext
7✔
719

720
    if file_extension == '.zip':
7✔
721
        with zipfile.ZipFile(file_path) as zf:
7✔
722
            members = zf.namelist()
7✔
723
            if len(members) != 1:
7✔
724
                raise RuntimeError('Cannot extract multiple files')
×
725
            o_name = hid + members[0]
7✔
726
            o_path = os.path.join(tmpdir, o_name)
7✔
727
            of_ext = _check_ext(o_path)
7✔
728
            if not os.path.exists(o_path):
7✔
729
                logger.info('Extracting {} to {}...'.format(bname, o_path))
7✔
730
                with open(o_path, 'wb') as f:
7✔
731
                    f.write(zf.read(members[0]))
7✔
732
    elif file_extension == '.gz':
4✔
733
        # Gzip files cannot be inspected. It's always only one file
734
        # Decide on its name
735
        o_name = hid + os.path.basename(filename)
3✔
736
        o_path = os.path.join(tmpdir, o_name)
3✔
737
        of_ext = _check_ext(o_path)
3✔
738
        if not os.path.exists(o_path):
3✔
739
            logger.info('Extracting {} to {}...'.format(bname, o_path))
3✔
740
            with gzip.GzipFile(file_path) as zf:
3✔
741
                with open(o_path, 'wb') as outfile:
3✔
742
                    for line in zf:
3✔
743
                        outfile.write(line)
3✔
744
    elif file_extension == '.bz2':
3✔
745
        # bzip2 files cannot be inspected. It's always only one file
746
        # Decide on its name
747
        o_name = hid + os.path.basename(filename)
3✔
748
        o_path = os.path.join(tmpdir, o_name)
3✔
749
        of_ext = _check_ext(o_path)
3✔
750
        if not os.path.exists(o_path):
3✔
751
            logger.info('Extracting {} to {}...'.format(bname, o_path))
3✔
752
            with bz2.open(file_path) as zf:
3✔
753
                with open(o_path, 'wb') as outfile:
3✔
754
                    for line in zf:
3✔
755
                        outfile.write(line)
3✔
756
    elif file_extension in ['.tar.gz', '.tar']:
1✔
757
        with tarfile.open(file_path) as zf:
1✔
758
            members = zf.getmembers()
1✔
759
            if len(members) != 1:
1✔
760
                raise RuntimeError('Cannot extract multiple files')
×
761
            o_name = hid + members[0].name
1✔
762
            o_path = os.path.join(tmpdir, o_name)
1✔
763
            of_ext = _check_ext(o_path)
1✔
764
            if not os.path.exists(o_path):
1✔
765
                logger.info('Extracting {} to {}...'.format(bname, o_path))
1✔
766
                with open(o_path, 'wb') as f:
1✔
767
                    f.write(zf.extractfile(members[0]).read())
1✔
768
    else:
769
        raise InvalidParamsError('Extension not recognized: '
1✔
770
                                 '{}'.format(file_extension))
771

772
    # Be sure we don't overfill the folder
773
    cfg.get_lru_handler(tmpdir, ending=of_ext).append(o_path)
7✔
774

775
    return o_path
7✔
776

777

778
def download_with_authentication(wwwfile, key):
9✔
779
    """ Uses credentials from a local .netrc file to download files
780

781
    This is function is currently used for TanDEM-X and ASTER
782

783
    Parameters
784
    ----------
785
    wwwfile : str
786
        path to the file to download
787
    key : str
788
        the machine to to look at in the .netrc file
789

790
    Returns
791
    -------
792

793
    """
794
    # Check the cache first. Use dummy download function to assure nothing is
795
    # tried to be downloaded without credentials:
796

797
    def _always_none(foo):
×
798
        return None
×
799

800
    cache_obj_name = _get_url_cache_name(wwwfile)
×
801
    dest_file = _verified_download_helper(cache_obj_name, _always_none)
×
802

803
    # Grab auth parameters
804
    if not dest_file:
×
805
        authfile = os.path.expanduser('~/.netrc')
×
806

807
        if not os.path.isfile(authfile):
×
808
            raise DownloadCredentialsMissingException(
×
809
                (authfile, ' does not exist. Add necessary credentials for ',
810
                 key, ' with `oggm_netrc_credentials. You may have to ',
811
                 'register at the respective service first.'))
812

813
        try:
×
814
            netrc(authfile).authenticators(key)[0]
×
815
        except TypeError:
×
816
            raise DownloadCredentialsMissingException(
×
817
                ('Credentials for ', key, ' are not in ', authfile, '. Add ',
818
                 'credentials for with `oggm_netrc_credentials`.'))
819

820
        dest_file = file_downloader(
×
821
            wwwfile, auth=(netrc(authfile).authenticators(key)[0],
822
                           netrc(authfile).authenticators(key)[2]))
823

824
    return dest_file
×
825

826

827
def download_oggm_files():
9✔
828
    with get_lock():
8✔
829
        return _download_oggm_files_unlocked()
8✔
830

831

832
def _download_oggm_files_unlocked():
9✔
833
    """Checks if the demo data is already on the cache and downloads it."""
834

835
    zip_url = 'https://github.com/%s/archive/%s.zip' % \
8✔
836
              (SAMPLE_DATA_GH_REPO, SAMPLE_DATA_COMMIT)
837
    odir = os.path.join(cfg.CACHE_DIR)
8✔
838
    sdir = os.path.join(cfg.CACHE_DIR,
8✔
839
                        'oggm-sample-data-%s' % SAMPLE_DATA_COMMIT)
840

841
    # download only if necessary
842
    if not os.path.exists(sdir):
8✔
843
        ofile = file_downloader(zip_url)
8✔
844
        with zipfile.ZipFile(ofile) as zf:
8✔
845
            zf.extractall(odir)
8✔
846
        assert os.path.isdir(sdir)
8✔
847

848
    # list of files for output
849
    out = dict()
8✔
850
    for root, directories, filenames in os.walk(sdir):
8✔
851
        for filename in filenames:
8✔
852
            if filename in out:
8✔
853
                # This was a stupid thing, and should not happen
854
                # TODO: duplicates in sample data...
855
                k = os.path.join(os.path.basename(root), filename)
8✔
856
                assert k not in out
8✔
857
                out[k] = os.path.join(root, filename)
8✔
858
            else:
859
                out[filename] = os.path.join(root, filename)
8✔
860

861
    return out
8✔
862

863

864
def _download_srtm_file(zone):
9✔
865
    with get_lock():
1✔
866
        return _download_srtm_file_unlocked(zone)
1✔
867

868

869
def _download_srtm_file_unlocked(zone):
9✔
870
    """Checks if the srtm data is in the directory and if not, download it.
871
    """
872

873
    # extract directory
874
    tmpdir = cfg.PATHS['tmp_dir']
1✔
875
    mkdir(tmpdir)
1✔
876
    outpath = os.path.join(tmpdir, 'srtm_' + zone + '.tif')
1✔
877

878
    # check if extracted file exists already
879
    if os.path.exists(outpath):
1✔
880
        return outpath
×
881

882
    # Did we download it yet?
883
    wwwfile = ('http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/'
1✔
884
               'TIFF/srtm_' + zone + '.zip')
885
    dest_file = file_downloader(wwwfile)
1✔
886

887
    # None means we tried hard but we couldn't find it
888
    if not dest_file:
1✔
889
        return None
×
890

891
    # ok we have to extract it
892
    if not os.path.exists(outpath):
1✔
893
        with zipfile.ZipFile(dest_file) as zf:
1✔
894
            zf.extractall(tmpdir)
1✔
895

896
    # See if we're good, don't overfill the tmp directory
897
    assert os.path.exists(outpath)
1✔
898
    cfg.get_lru_handler(tmpdir).append(outpath)
1✔
899
    return outpath
1✔
900

901

902
def _download_nasadem_file(zone):
9✔
903
    with get_lock():
1✔
904
        return _download_nasadem_file_unlocked(zone)
1✔
905

906

907
def _download_nasadem_file_unlocked(zone):
9✔
908
    """Checks if the NASADEM data is in the directory and if not, download it.
909
    """
910

911
    # extract directory
912
    tmpdir = cfg.PATHS['tmp_dir']
1✔
913
    mkdir(tmpdir)
1✔
914
    wwwfile = ('https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/'
1✔
915
               '2000.02.11/NASADEM_HGT_{}.zip'.format(zone))
916
    demfile = '{}.hgt'.format(zone)
1✔
917
    outpath = os.path.join(tmpdir, demfile)
1✔
918

919
    # check if extracted file exists already
920
    if os.path.exists(outpath):
1✔
921
        return outpath
×
922

923
    # Did we download it yet?
924
    dest_file = file_downloader(wwwfile)
1✔
925

926
    # None means we tried hard but we couldn't find it
927
    if not dest_file:
1✔
928
        return None
×
929

930
    # ok we have to extract it
931
    if not os.path.exists(outpath):
1✔
932
        with zipfile.ZipFile(dest_file) as zf:
1✔
933
            zf.extract(demfile, path=tmpdir)
1✔
934

935
    # See if we're good, don't overfill the tmp directory
936
    assert os.path.exists(outpath)
1✔
937
    cfg.get_lru_handler(tmpdir).append(outpath)
1✔
938
    return outpath
1✔
939

940

941
def _download_tandem_file(zone):
9✔
942
    with get_lock():
1✔
943
        return _download_tandem_file_unlocked(zone)
1✔
944

945

946
def _download_tandem_file_unlocked(zone):
9✔
947
    """Checks if the tandem data is in the directory and if not, download it.
948
    """
949

950
    # extract directory
951
    tmpdir = cfg.PATHS['tmp_dir']
1✔
952
    mkdir(tmpdir)
1✔
953
    bname = zone.split('/')[-1] + '_DEM.tif'
1✔
954
    wwwfile = ('https://download.geoservice.dlr.de/TDM90/files/DEM/'
1✔
955
               '{}.zip'.format(zone))
956
    outpath = os.path.join(tmpdir, bname)
1✔
957

958
    # check if extracted file exists already
959
    if os.path.exists(outpath):
1✔
960
        return outpath
×
961

962
    dest_file = download_with_authentication(wwwfile, 'geoservice.dlr.de')
1✔
963

964
    # That means we tried hard but we couldn't find it
965
    if not dest_file:
1✔
966
        return None
×
967
    elif not zipfile.is_zipfile(dest_file):
1✔
968
        # If the TanDEM-X tile does not exist, a invalid file is created.
969
        # See https://github.com/OGGM/oggm/issues/893 for more details
970
        return None
1✔
971

972
    # ok we have to extract it
973
    if not os.path.exists(outpath):
1✔
974
        with zipfile.ZipFile(dest_file) as zf:
1✔
975
            for fn in zf.namelist():
1✔
976
                if 'DEM/' + bname in fn:
1✔
977
                    break
×
978
            with open(outpath, 'wb') as fo:
1✔
979
                fo.write(zf.read(fn))
1✔
980

981
    # See if we're good, don't overfill the tmp directory
982
    assert os.path.exists(outpath)
1✔
983
    cfg.get_lru_handler(tmpdir).append(outpath)
1✔
984
    return outpath
1✔
985

986

987
def _download_dem3_viewpano(zone):
9✔
988
    with get_lock():
1✔
989
        return _download_dem3_viewpano_unlocked(zone)
1✔
990

991

992
def _download_dem3_viewpano_unlocked(zone):
9✔
993
    """Checks if the DEM3 data is in the directory and if not, download it.
994
    """
995

996
    # extract directory
997
    tmpdir = cfg.PATHS['tmp_dir']
1✔
998
    mkdir(tmpdir)
1✔
999
    outpath = os.path.join(tmpdir, zone + '.tif')
1✔
1000
    extract_dir = os.path.join(tmpdir, 'tmp_' + zone)
1✔
1001
    mkdir(extract_dir, reset=True)
1✔
1002

1003
    # check if extracted file exists already
1004
    if os.path.exists(outpath):
1✔
1005
        return outpath
×
1006

1007
    # OK, so see if downloaded already
1008
    # some files have a newer version 'v2'
1009
    if zone in ['R33', 'R34', 'R35', 'R36', 'R37', 'R38', 'Q32', 'Q33', 'Q34',
1✔
1010
                'Q35', 'Q36', 'Q37', 'Q38', 'Q39', 'Q40', 'P31', 'P32', 'P33',
1011
                'P34', 'P35', 'P36', 'P37', 'P38', 'P39', 'P40']:
1012
        ifile = 'http://viewfinderpanoramas.org/dem3/' + zone + 'v2.zip'
×
1013
    elif zone in DEM3REG.keys():
1✔
1014
        # We prepared these files as tif already
1015
        ifile = ('https://cluster.klima.uni-bremen.de/~oggm/dem/'
1✔
1016
                 'DEM3_MERGED/{}.tif'.format(zone))
1017
        return file_downloader(ifile)
1✔
1018
    else:
1019
        ifile = 'http://viewfinderpanoramas.org/dem3/' + zone + '.zip'
1✔
1020

1021
    dfile = file_downloader(ifile)
1✔
1022

1023
    # None means we tried hard but we couldn't find it
1024
    if not dfile:
1✔
1025
        return None
×
1026

1027
    # ok we have to extract it
1028
    with zipfile.ZipFile(dfile) as zf:
1✔
1029
        zf.extractall(extract_dir)
1✔
1030

1031
    # Serious issue: sometimes, if a southern hemisphere URL is queried for
1032
    # download and there is none, a NH zip file is downloaded.
1033
    # Example: http://viewfinderpanoramas.org/dem3/SN29.zip yields N29!
1034
    # BUT: There are southern hemisphere files that download properly. However,
1035
    # the unzipped folder has the file name of
1036
    # the northern hemisphere file. Some checks if correct file exists:
1037
    if len(zone) == 4 and zone.startswith('S'):
1✔
1038
        zonedir = os.path.join(extract_dir, zone[1:])
×
1039
    else:
1040
        zonedir = os.path.join(extract_dir, zone)
1✔
1041
    globlist = glob.glob(os.path.join(zonedir, '*.hgt'))
1✔
1042

1043
    # take care of the special file naming cases
1044
    if zone in DEM3REG.keys():
1✔
1045
        globlist = glob.glob(os.path.join(extract_dir, '*', '*.hgt'))
×
1046

1047
    if not globlist:
1✔
1048
        # Final resort
1049
        globlist = (findfiles(extract_dir, '.hgt') or
×
1050
                    findfiles(extract_dir, '.HGT'))
1051
        if not globlist:
×
1052
            raise RuntimeError("We should have some files here, but we don't")
×
1053

1054
    # merge the single HGT files (can be a bit ineffective, because not every
1055
    # single file might be exactly within extent...)
1056
    rfiles = [rasterio.open(s) for s in globlist]
1✔
1057
    dest, output_transform = merge_tool(rfiles)
1✔
1058
    profile = rfiles[0].profile
1✔
1059
    if 'affine' in profile:
1✔
1060
        profile.pop('affine')
×
1061
    profile['transform'] = output_transform
1✔
1062
    profile['height'] = dest.shape[1]
1✔
1063
    profile['width'] = dest.shape[2]
1✔
1064
    profile['driver'] = 'GTiff'
1✔
1065
    with rasterio.open(outpath, 'w', **profile) as dst:
1✔
1066
        dst.write(dest)
1✔
1067
    for rf in rfiles:
1✔
1068
        rf.close()
1✔
1069

1070
    # delete original files to spare disk space
1071
    for s in globlist:
1✔
1072
        os.remove(s)
1✔
1073
    del_empty_dirs(tmpdir)
1✔
1074

1075
    # See if we're good, don't overfill the tmp directory
1076
    assert os.path.exists(outpath)
1✔
1077
    cfg.get_lru_handler(tmpdir).append(outpath)
1✔
1078
    return outpath
1✔
1079

1080

1081
def _download_aster_file(zone):
9✔
1082
    with get_lock():
1✔
1083
        return _download_aster_file_unlocked(zone)
1✔
1084

1085

1086
def _download_aster_file_unlocked(zone):
9✔
1087
    """Checks if the ASTER data is in the directory and if not, download it.
1088
    """
1089

1090
    # extract directory
1091
    tmpdir = cfg.PATHS['tmp_dir']
1✔
1092
    mkdir(tmpdir)
1✔
1093
    wwwfile = ('https://e4ftl01.cr.usgs.gov/ASTT/ASTGTM.003/'
1✔
1094
               '2000.03.01/{}.zip'.format(zone))
1095
    outpath = os.path.join(tmpdir, zone + '_dem.tif')
1✔
1096

1097
    # check if extracted file exists already
1098
    if os.path.exists(outpath):
1✔
1099
        return outpath
×
1100

1101
    # download from NASA Earthdata with credentials
1102
    dest_file = download_with_authentication(wwwfile, 'urs.earthdata.nasa.gov')
1✔
1103

1104
    # That means we tried hard but we couldn't find it
1105
    if not dest_file:
1✔
1106
        return None
×
1107

1108
    # ok we have to extract it
1109
    if not os.path.exists(outpath):
1✔
1110
        with zipfile.ZipFile(dest_file) as zf:
1✔
1111
            zf.extractall(tmpdir)
1✔
1112

1113
    # See if we're good, don't overfill the tmp directory
1114
    assert os.path.exists(outpath)
1✔
1115
    cfg.get_lru_handler(tmpdir).append(outpath)
1✔
1116
    return outpath
1✔
1117

1118

1119
def _download_topo_file_from_cluster(fname):
9✔
1120
    with get_lock():
×
1121
        return _download_topo_file_from_cluster_unlocked(fname)
×
1122

1123

1124
def _download_topo_file_from_cluster_unlocked(fname):
9✔
1125
    """Checks if the special topo data is in the directory and if not,
1126
    download it from the cluster.
1127
    """
1128

1129
    # extract directory
1130
    tmpdir = cfg.PATHS['tmp_dir']
×
1131
    mkdir(tmpdir)
×
1132
    outpath = os.path.join(tmpdir, fname)
×
1133

1134
    url = 'https://cluster.klima.uni-bremen.de/data/dems/'
×
1135
    url += fname + '.zip'
×
1136
    dfile = file_downloader(url)
×
1137

1138
    if not os.path.exists(outpath):
×
1139
        logger.info('Extracting ' + fname + '.zip to ' + outpath + '...')
×
1140
        with zipfile.ZipFile(dfile) as zf:
×
1141
            zf.extractall(tmpdir)
×
1142

1143
    # See if we're good, don't overfill the tmp directory
1144
    assert os.path.exists(outpath)
×
1145
    cfg.get_lru_handler(tmpdir).append(outpath)
×
1146
    return outpath
×
1147

1148

1149
def _download_copdem_file(cppfile, tilename, source):
9✔
1150
    with get_lock():
1✔
1151
        return _download_copdem_file_unlocked(cppfile, tilename, source)
1✔
1152

1153

1154
def _download_copdem_file_unlocked(cppfile, tilename, source):
9✔
1155
    """Checks if Copernicus DEM file is in the directory, if not download it.
1156

1157
    cppfile : name of the tarfile to download
1158
    tilename : name of folder and tif file within the cppfile
1159
    source : either 'COPDEM90' or 'COPDEM30'
1160

1161
    """
1162

1163
    # extract directory
1164
    tmpdir = cfg.PATHS['tmp_dir']
1✔
1165
    mkdir(tmpdir)
1✔
1166

1167
    # tarfiles are extracted in directories per each tile
1168
    fpath = '{0}_DEM.tif'.format(tilename)
1✔
1169
    demfile = os.path.join(tmpdir, fpath)
1✔
1170

1171
    # check if extracted file exists already
1172
    if os.path.exists(demfile):
1✔
1173
        return demfile
×
1174

1175
    # Did we download it yet?
1176
    ftpfile = ('ftps://cdsdata.copernicus.eu:990/' +
1✔
1177
               'datasets/COP-DEM_GLO-{}-DGED/2021_1/'.format(source[-2:]) +
1178
               cppfile)
1179

1180
    dest_file = download_with_authentication(ftpfile,
1✔
1181
                                             'spacedata.copernicus.eu')
1182

1183
    # None means we tried hard but we couldn't find it
1184
    if not dest_file:
1✔
1185
        return None
×
1186

1187
    # ok we have to extract it
1188
    if not os.path.exists(demfile):
1✔
1189
        tiffile = os.path.join(tilename, 'DEM', fpath)
1✔
1190
        with tarfile.open(dest_file) as tf:
1✔
1191
            tmember = tf.getmember(tiffile)
1✔
1192
            # do not extract the full path of the file
1193
            tmember.name = os.path.basename(tf.getmember(tiffile).name)
1✔
1194
            tf.extract(tmember, tmpdir)
1✔
1195

1196
    # See if we're good, don't overfill the tmp directory
1197
    assert os.path.exists(demfile)
1✔
1198
    cfg.get_lru_handler(tmpdir).append(demfile)
1✔
1199

1200
    return demfile
1✔
1201

1202

1203
def _download_aw3d30_file(zone):
9✔
1204
    with get_lock():
1✔
1205
        return _download_aw3d30_file_unlocked(zone)
1✔
1206

1207

1208
def _download_aw3d30_file_unlocked(fullzone):
9✔
1209
    """Checks if the AW3D30 data is in the directory and if not, download it.
1210
    """
1211

1212
    # extract directory
1213
    tmpdir = cfg.PATHS['tmp_dir']
1✔
1214
    mkdir(tmpdir)
1✔
1215

1216
    # tarfiles are extracted in directories per each tile
1217
    tile = fullzone.split('/')[1]
1✔
1218
    demfile = os.path.join(tmpdir, tile, tile + '_AVE_DSM.tif')
1✔
1219

1220
    # check if extracted file exists already
1221
    if os.path.exists(demfile):
1✔
1222
        return demfile
×
1223

1224
    # Did we download it yet?
1225
    ftpfile = ('ftp://ftp.eorc.jaxa.jp/pub/ALOS/ext1/AW3D30/release_v1804/'
1✔
1226
               + fullzone + '.tar.gz')
1227
    try:
1✔
1228
        dest_file = file_downloader(ftpfile, timeout=180)
1✔
1229
    except urllib.error.URLError:
×
1230
        # This error is raised if file is not available, could be water
1231
        return None
×
1232

1233
    # None means we tried hard but we couldn't find it
1234
    if not dest_file:
1✔
1235
        return None
×
1236

1237
    # ok we have to extract it
1238
    if not os.path.exists(demfile):
1✔
1239
        from oggm.utils import robust_tar_extract
1✔
1240
        dempath = os.path.dirname(demfile)
1✔
1241
        robust_tar_extract(dest_file, dempath)
1✔
1242

1243
    # See if we're good, don't overfill the tmp directory
1244
    assert os.path.exists(demfile)
1✔
1245
    # this tarfile contains several files
1246
    for file in os.listdir(dempath):
1✔
1247
        cfg.get_lru_handler(tmpdir).append(os.path.join(dempath, file))
1✔
1248
    return demfile
1✔
1249

1250

1251
def _download_mapzen_file(zone):
9✔
1252
    with get_lock():
1✔
1253
        return _download_mapzen_file_unlocked(zone)
1✔
1254

1255

1256
def _download_mapzen_file_unlocked(zone):
9✔
1257
    """Checks if the mapzen data is in the directory and if not, download it.
1258
    """
1259
    bucket = 'elevation-tiles-prod'
1✔
1260
    prefix = 'geotiff'
1✔
1261
    url = 'http://s3.amazonaws.com/%s/%s/%s' % (bucket, prefix, zone)
1✔
1262

1263
    # That's all
1264
    return file_downloader(url, timeout=180)
1✔
1265

1266

1267
def get_prepro_gdir(rgi_version, rgi_id, border, prepro_level, base_url=None):
9✔
1268
    with get_lock():
2✔
1269
        return _get_prepro_gdir_unlocked(rgi_version, rgi_id, border,
2✔
1270
                                         prepro_level, base_url=base_url)
1271

1272

1273
def get_prepro_base_url(base_url=None, rgi_version=None, border=None,
9✔
1274
                        prepro_level=None):
1275
    """Extended base url where to find the desired gdirs."""
1276

1277
    if base_url is None:
2✔
1278
        raise InvalidParamsError('Starting with v1.6, users now have to '
1✔
1279
                                 'explicitly indicate the url they want '
1280
                                 'to start from.')
1281

1282
    if not base_url.endswith('/'):
2✔
1283
        base_url += '/'
×
1284

1285
    if rgi_version is None:
2✔
1286
        rgi_version = cfg.PARAMS['rgi_version']
1✔
1287

1288
    if border is None:
2✔
1289
        border = cfg.PARAMS['border']
1✔
1290

1291
    url = base_url
2✔
1292
    url += 'RGI{}/'.format(rgi_version)
2✔
1293
    url += 'b_{:03d}/'.format(int(border))
2✔
1294
    url += 'L{:d}/'.format(prepro_level)
2✔
1295
    return url
2✔
1296

1297

1298
def _get_prepro_gdir_unlocked(rgi_version, rgi_id, border, prepro_level,
9✔
1299
                              base_url=None):
1300

1301
    url = get_prepro_base_url(rgi_version=rgi_version, border=border,
2✔
1302
                              prepro_level=prepro_level, base_url=base_url)
1303
    url += '{}/{}.tar' .format(rgi_id[:8], rgi_id[:11])
2✔
1304
    tar_base = file_downloader(url)
2✔
1305
    if tar_base is None:
2✔
1306
        raise RuntimeError('Could not find file at ' + url)
×
1307

1308
    return tar_base
2✔
1309

1310

1311
def get_geodetic_mb_dataframe(file_path=None):
9✔
1312
    """Fetches the reference geodetic dataframe for calibration.
1313

1314
    Currently that's the data from Hughonnet et al 2021, corrected for
1315
    outliers and with void filled. The data preparation script is
1316
    available at
1317
    https://nbviewer.jupyter.org/urls/cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb/convert.ipynb
1318

1319
    Parameters
1320
    ----------
1321
    file_path : str
1322
        in case you have your own file to parse (check the format first!)
1323

1324
    Returns
1325
    -------
1326
    a DataFrame with the data.
1327
    """
1328

1329
    # fetch the file online or read custom file
1330
    if file_path is None:
3✔
1331
        base_url = 'https://cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb/'
3✔
1332
        file_name = 'hugonnet_2021_ds_rgi60_pergla_rates_10_20_worldwide_filled.hdf'
3✔
1333
        file_path = file_downloader(base_url + file_name)
3✔
1334

1335
    # Did we open it yet?
1336
    if file_path in cfg.DATA:
3✔
1337
        return cfg.DATA[file_path]
2✔
1338

1339
    # If not let's go
1340
    extension = os.path.splitext(file_path)[1]
3✔
1341
    if extension == '.csv':
3✔
1342
        df = pd.read_csv(file_path, index_col=0)
×
1343
    elif extension == '.hdf':
3✔
1344
        df = pd.read_hdf(file_path)
3✔
1345

1346
    # Check for missing data (old files)
1347
    if len(df.loc[df['dmdtda'].isnull()]) > 0:
3✔
1348
        raise InvalidParamsError('The reference file you are using has missing '
×
1349
                                 'data and is probably outdated (sorry for '
1350
                                 'that). Delete the file at '
1351
                                 f'{file_path} and start again.')
1352
    cfg.DATA[file_path] = df
3✔
1353
    return df
3✔
1354

1355

1356
def srtm_zone(lon_ex, lat_ex):
9✔
1357
    """Returns a list of SRTM zones covering the desired extent.
1358
    """
1359

1360
    # SRTM are sorted in tiles of 5 degrees
1361
    srtm_x0 = -180.
1✔
1362
    srtm_y0 = 60.
1✔
1363
    srtm_dx = 5.
1✔
1364
    srtm_dy = -5.
1✔
1365

1366
    # quick n dirty solution to be sure that we will cover the whole range
1367
    mi, ma = np.min(lon_ex), np.max(lon_ex)
1✔
1368
    # int() to avoid Deprec warning:
1369
    lon_ex = np.linspace(mi, ma, int(np.ceil((ma - mi) + 3)))
1✔
1370
    mi, ma = np.min(lat_ex), np.max(lat_ex)
1✔
1371
    # int() to avoid Deprec warning
1372
    lat_ex = np.linspace(mi, ma, int(np.ceil((ma - mi) + 3)))
1✔
1373

1374
    zones = []
1✔
1375
    for lon in lon_ex:
1✔
1376
        for lat in lat_ex:
1✔
1377
            dx = lon - srtm_x0
1✔
1378
            dy = lat - srtm_y0
1✔
1379
            if dy > 0:
1✔
1380
                continue
×
1381
            zx = np.ceil(dx / srtm_dx)
1✔
1382
            zy = np.ceil(dy / srtm_dy)
1✔
1383
            zones.append('{:02.0f}_{:02.0f}'.format(zx, zy))
1✔
1384
    return list(sorted(set(zones)))
1✔
1385

1386

1387
def _tandem_path(lon_tile, lat_tile):
9✔
1388

1389
    # OK we have a proper tile now
1390
    # This changed in December 2022
1391

1392
    # First folder level is sorted from S to N
1393
    level_0 = 'S' if lat_tile < 0 else 'N'
1✔
1394
    level_0 += '{:02d}'.format(abs(lat_tile))
1✔
1395

1396
    # Second folder level is sorted from W to E, but in 10 steps
1397
    level_1 = 'W' if lon_tile < 0 else 'E'
1✔
1398
    level_1 += '{:03d}'.format(divmod(abs(lon_tile), 10)[0] * 10)
1✔
1399

1400
    # Level 2 is formatting, but depends on lat
1401
    level_2 = 'W' if lon_tile < 0 else 'E'
1✔
1402
    if abs(lat_tile) <= 60:
1✔
1403
        level_2 += '{:03d}'.format(abs(lon_tile))
1✔
1404
    elif abs(lat_tile) <= 80:
1✔
1405
        level_2 += '{:03d}'.format(divmod(abs(lon_tile), 2)[0] * 2)
1✔
1406
    else:
1407
        level_2 += '{:03d}'.format(divmod(abs(lon_tile), 4)[0] * 4)
1✔
1408

1409
    # Final path
1410
    out = (level_0 + '/' + level_1 + '/' +
1✔
1411
           'TDM1_DEM__30_{}{}'.format(level_0, level_2))
1412
    return out
1✔
1413

1414

1415
def tandem_zone(lon_ex, lat_ex):
9✔
1416
    """Returns a list of TanDEM-X zones covering the desired extent.
1417
    """
1418

1419
    # Files are one by one tiles, so lets loop over them
1420
    # For higher lats they are stored in steps of 2 and 4. My code below
1421
    # is probably giving more files than needed but better safe than sorry
1422
    lat_tiles = np.arange(np.floor(lat_ex[0]), np.ceil(lat_ex[1]+1e-9),
1✔
1423
                          dtype=int)
1424
    zones = []
1✔
1425
    for lat in lat_tiles:
1✔
1426
        if abs(lat) < 60:
1✔
1427
            l0 = np.floor(lon_ex[0])
1✔
1428
            l1 = np.floor(lon_ex[1])
1✔
1429
        elif abs(lat) < 80:
1✔
1430
            l0 = divmod(lon_ex[0], 2)[0] * 2
1✔
1431
            l1 = divmod(lon_ex[1], 2)[0] * 2
1✔
1432
        elif abs(lat) < 90:
1✔
1433
            l0 = divmod(lon_ex[0], 4)[0] * 4
1✔
1434
            l1 = divmod(lon_ex[1], 4)[0] * 4
1✔
1435
        lon_tiles = np.arange(l0, l1+1, dtype=int)
1✔
1436
        for lon in lon_tiles:
1✔
1437
            zones.append(_tandem_path(lon, lat))
1✔
1438
    return list(sorted(set(zones)))
1✔
1439

1440

1441
def _aw3d30_path(lon_tile, lat_tile):
9✔
1442

1443
    # OK we have a proper tile now
1444

1445
    # Folders are sorted with N E S W in 5 degree steps
1446
    # But in N and E the lower boundary is indicated
1447
    # e.g. N060 contains N060 - N064
1448
    # e.g. E000 contains E000 - E004
1449
    # but S and W indicate the upper boundary:
1450
    # e.g. S010 contains S006 - S010
1451
    # e.g. W095 contains W091 - W095
1452

1453
    # get letters
1454
    ns = 'S' if lat_tile < 0 else 'N'
1✔
1455
    ew = 'W' if lon_tile < 0 else 'E'
1✔
1456

1457
    # get lat/lon
1458
    lon = abs(5 * np.floor(lon_tile/5))
1✔
1459
    lat = abs(5 * np.floor(lat_tile/5))
1✔
1460

1461
    folder = '%s%.3d%s%.3d' % (ns, lat, ew, lon)
1✔
1462
    filename = '%s%.3d%s%.3d' % (ns, abs(lat_tile), ew, abs(lon_tile))
1✔
1463

1464
    # Final path
1465
    out = folder + '/' + filename
1✔
1466
    return out
1✔
1467

1468

1469
def aw3d30_zone(lon_ex, lat_ex):
9✔
1470
    """Returns a list of AW3D30 zones covering the desired extent.
1471
    """
1472

1473
    # Files are one by one tiles, so lets loop over them
1474
    lon_tiles = np.arange(np.floor(lon_ex[0]), np.ceil(lon_ex[1]+1e-9),
1✔
1475
                          dtype=int)
1476
    lat_tiles = np.arange(np.floor(lat_ex[0]), np.ceil(lat_ex[1]+1e-9),
1✔
1477
                          dtype=int)
1478
    zones = []
1✔
1479
    for lon in lon_tiles:
1✔
1480
        for lat in lat_tiles:
1✔
1481
            zones.append(_aw3d30_path(lon, lat))
1✔
1482
    return list(sorted(set(zones)))
1✔
1483

1484

1485
def _extent_to_polygon(lon_ex, lat_ex, to_crs=None):
9✔
1486

1487
    if lon_ex[0] == lon_ex[1] and lat_ex[0] == lat_ex[1]:
1✔
1488
        out = shpg.Point(lon_ex[0], lat_ex[0])
1✔
1489
    else:
1490
        x = [lon_ex[0], lon_ex[1], lon_ex[1], lon_ex[0], lon_ex[0]]
1✔
1491
        y = [lat_ex[0], lat_ex[0], lat_ex[1], lat_ex[1], lat_ex[0]]
1✔
1492
        out = shpg.Polygon(np.array((x, y)).T)
1✔
1493
    if to_crs is not None:
1✔
1494
        out = salem.transform_geometry(out, to_crs=to_crs)
1✔
1495
    return out
1✔
1496

1497

1498
def arcticdem_zone(lon_ex, lat_ex):
9✔
1499
    """Returns a list of Arctic-DEM zones covering the desired extent.
1500
    """
1501

1502
    gdf = gpd.read_file(get_demo_file('ArcticDEM_Tile_Index_Rel7_by_tile.shp'))
1✔
1503
    p = _extent_to_polygon(lon_ex, lat_ex, to_crs=gdf.crs)
1✔
1504
    gdf = gdf.loc[gdf.intersects(p)]
1✔
1505
    return gdf.tile.values if len(gdf) > 0 else []
1✔
1506

1507

1508
def rema_zone(lon_ex, lat_ex):
9✔
1509
    """Returns a list of REMA-DEM zones covering the desired extent.
1510
    """
1511

1512
    gdf = gpd.read_file(get_demo_file('REMA_Tile_Index_Rel1.1.shp'))
1✔
1513
    p = _extent_to_polygon(lon_ex, lat_ex, to_crs=gdf.crs)
1✔
1514
    gdf = gdf.loc[gdf.intersects(p)]
1✔
1515
    return gdf.tile.values if len(gdf) > 0 else []
1✔
1516

1517

1518
def alaska_dem_zone(lon_ex, lat_ex):
9✔
1519
    """Returns a list of Alaska-DEM zones covering the desired extent.
1520
    """
1521

1522
    gdf = gpd.read_file(get_demo_file('Alaska_albers_V3_tiles.shp'))
1✔
1523
    p = _extent_to_polygon(lon_ex, lat_ex, to_crs=gdf.crs)
1✔
1524
    gdf = gdf.loc[gdf.intersects(p)]
1✔
1525
    return gdf.tile.values if len(gdf) > 0 else []
1✔
1526

1527

1528
def copdem_zone(lon_ex, lat_ex, source):
9✔
1529
    """Returns a list of Copernicus DEM tarfile and tilename tuples
1530
    """
1531

1532
    # because we use both meters and arc secs in our filenames...
1533
    if source[-2:] == '90':
1✔
1534
        asec = '30'
1✔
1535
    elif source[-2:] == '30':
1✔
1536
        asec = '10'
1✔
1537
    else:
1538
        raise InvalidDEMError('COPDEM Version not valid.')
×
1539

1540
    # either reuse or load lookup table
1541
    if source in cfg.DATA:
1✔
1542
        df = cfg.DATA[source]
1✔
1543
    else:
1544
        df = pd.read_csv(get_demo_file('{}_2021_1.csv'.format(source.lower())))
1✔
1545
        cfg.DATA[source] = df
1✔
1546

1547
    # adding small buffer for unlikely case where one lon/lat_ex == xx.0
1548
    lons = np.arange(np.floor(lon_ex[0]-1e-9), np.ceil(lon_ex[1]+1e-9))
1✔
1549
    lats = np.arange(np.floor(lat_ex[0]-1e-9), np.ceil(lat_ex[1]+1e-9))
1✔
1550

1551
    flist = []
1✔
1552
    for lat in lats:
1✔
1553
        # north or south?
1554
        ns = 'S' if lat < 0 else 'N'
1✔
1555
        for lon in lons:
1✔
1556
            # east or west?
1557
            ew = 'W' if lon < 0 else 'E'
1✔
1558
            lat_str = '{}{:02.0f}'.format(ns, abs(lat))
1✔
1559
            lon_str = '{}{:03.0f}'.format(ew, abs(lon))
1✔
1560
            try:
1✔
1561
                filename = df.loc[(df['Long'] == lon_str) &
1✔
1562
                                  (df['Lat'] == lat_str)]['CPP filename'].iloc[0]
1563
                flist.append((filename,
1✔
1564
                              'Copernicus_DSM_{}_{}_00_{}_00'.format(asec, lat_str, lon_str)))
1565
            except IndexError:
×
1566
                # COPDEM is global, if we miss tiles it is probably in the ocean
1567
                pass
×
1568
    return flist
1✔
1569

1570

1571
def dem3_viewpano_zone(lon_ex, lat_ex):
9✔
1572
    """Returns a list of DEM3 zones covering the desired extent.
1573

1574
    http://viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm
1575
    """
1576

1577
    for _f in DEM3REG.keys():
1✔
1578

1579
        if (np.min(lon_ex) >= DEM3REG[_f][0]) and \
1✔
1580
           (np.max(lon_ex) <= DEM3REG[_f][1]) and \
1581
           (np.min(lat_ex) >= DEM3REG[_f][2]) and \
1582
           (np.max(lat_ex) <= DEM3REG[_f][3]):
1583

1584
            # test some weird inset files in Antarctica
1585
            if (np.min(lon_ex) >= -91.) and (np.max(lon_ex) <= -90.) and \
1✔
1586
               (np.min(lat_ex) >= -72.) and (np.max(lat_ex) <= -68.):
1587
                return ['SR15']
×
1588

1589
            elif (np.min(lon_ex) >= -47.) and (np.max(lon_ex) <= -43.) and \
1✔
1590
                 (np.min(lat_ex) >= -61.) and (np.max(lat_ex) <= -60.):
1591
                return ['SP23']
×
1592

1593
            elif (np.min(lon_ex) >= 162.) and (np.max(lon_ex) <= 165.) and \
1✔
1594
                 (np.min(lat_ex) >= -68.) and (np.max(lat_ex) <= -66.):
1595
                return ['SQ58']
×
1596

1597
            # test some rogue Greenland tiles as well
1598
            elif (np.min(lon_ex) >= -72.) and (np.max(lon_ex) <= -66.) and \
1✔
1599
                 (np.min(lat_ex) >= 76.) and (np.max(lat_ex) <= 80.):
1600
                return ['T19']
×
1601

1602
            elif (np.min(lon_ex) >= -72.) and (np.max(lon_ex) <= -66.) and \
1✔
1603
                 (np.min(lat_ex) >= 80.) and (np.max(lat_ex) <= 83.):
1604
                return ['U19']
×
1605

1606
            elif (np.min(lon_ex) >= -66.) and (np.max(lon_ex) <= -60.) and \
1✔
1607
                 (np.min(lat_ex) >= 80.) and (np.max(lat_ex) <= 83.):
1608
                return ['U20']
×
1609

1610
            elif (np.min(lon_ex) >= -60.) and (np.max(lon_ex) <= -54.) and \
1✔
1611
                 (np.min(lat_ex) >= 80.) and (np.max(lat_ex) <= 83.):
1612
                return ['U21']
×
1613

1614
            elif (np.min(lon_ex) >= -54.) and (np.max(lon_ex) <= -48.) and \
1✔
1615
                 (np.min(lat_ex) >= 80.) and (np.max(lat_ex) <= 83.):
1616
                return ['U22']
×
1617

1618
            elif (np.min(lon_ex) >= -25.) and (np.max(lon_ex) <= -13.) and \
1✔
1619
                 (np.min(lat_ex) >= 63.) and (np.max(lat_ex) <= 67.):
1620
                return ['ISL']
1✔
1621

1622
            else:
1623
                return [_f]
1✔
1624

1625
    # if the tile doesn't have a special name, its name can be found like this:
1626
    # corrected SRTMs are sorted in tiles of 6 deg longitude and 4 deg latitude
1627
    srtm_x0 = -180.
1✔
1628
    srtm_y0 = 0.
1✔
1629
    srtm_dx = 6.
1✔
1630
    srtm_dy = 4.
1✔
1631

1632
    # quick n dirty solution to be sure that we will cover the whole range
1633
    mi, ma = np.min(lon_ex), np.max(lon_ex)
1✔
1634
    # TODO: Fabien, find out what Johannes wanted with this +3
1635
    # +3 is just for the number to become still a bit larger
1636
    # int() to avoid Deprec warning
1637
    lon_ex = np.linspace(mi, ma, int(np.ceil((ma - mi) / srtm_dy) + 3))
1✔
1638
    mi, ma = np.min(lat_ex), np.max(lat_ex)
1✔
1639
    # int() to avoid Deprec warning
1640
    lat_ex = np.linspace(mi, ma, int(np.ceil((ma - mi) / srtm_dx) + 3))
1✔
1641

1642
    zones = []
1✔
1643
    for lon in lon_ex:
1✔
1644
        for lat in lat_ex:
1✔
1645
            dx = lon - srtm_x0
1✔
1646
            dy = lat - srtm_y0
1✔
1647
            zx = np.ceil(dx / srtm_dx)
1✔
1648
            # convert number to letter
1649
            zy = chr(int(abs(dy / srtm_dy)) + ord('A'))
1✔
1650
            if lat >= 0:
1✔
1651
                zones.append('%s%02.0f' % (zy, zx))
1✔
1652
            else:
1653
                zones.append('S%s%02.0f' % (zy, zx))
×
1654
    return list(sorted(set(zones)))
1✔
1655

1656

1657
def aster_zone(lon_ex, lat_ex):
9✔
1658
    """Returns a list of ASTGTMV3 zones covering the desired extent.
1659

1660
    ASTER v3 tiles are 1 degree x 1 degree
1661
    N50 contains 50 to 50.9
1662
    E10 contains 10 to 10.9
1663
    S70 contains -69.99 to -69.0
1664
    W20 contains -19.99 to -19.0
1665
    """
1666

1667
    # adding small buffer for unlikely case where one lon/lat_ex == xx.0
1668
    lons = np.arange(np.floor(lon_ex[0]-1e-9), np.ceil(lon_ex[1]+1e-9))
1✔
1669
    lats = np.arange(np.floor(lat_ex[0]-1e-9), np.ceil(lat_ex[1]+1e-9))
1✔
1670

1671
    zones = []
1✔
1672
    for lat in lats:
1✔
1673
        # north or south?
1674
        ns = 'S' if lat < 0 else 'N'
1✔
1675
        for lon in lons:
1✔
1676
            # east or west?
1677
            ew = 'W' if lon < 0 else 'E'
1✔
1678
            filename = 'ASTGTMV003_{}{:02.0f}{}{:03.0f}'.format(ns, abs(lat),
1✔
1679
                                                                ew, abs(lon))
1680
            zones.append(filename)
1✔
1681
    return list(sorted(set(zones)))
1✔
1682

1683

1684
def nasadem_zone(lon_ex, lat_ex):
9✔
1685
    """Returns a list of NASADEM zones covering the desired extent.
1686

1687
    NASADEM tiles are 1 degree x 1 degree
1688
    N50 contains 50 to 50.9
1689
    E10 contains 10 to 10.9
1690
    S70 contains -69.99 to -69.0
1691
    W20 contains -19.99 to -19.0
1692
    """
1693

1694
    # adding small buffer for unlikely case where one lon/lat_ex == xx.0
1695
    lons = np.arange(np.floor(lon_ex[0]-1e-9), np.ceil(lon_ex[1]+1e-9))
1✔
1696
    lats = np.arange(np.floor(lat_ex[0]-1e-9), np.ceil(lat_ex[1]+1e-9))
1✔
1697

1698
    zones = []
1✔
1699
    for lat in lats:
1✔
1700
        # north or south?
1701
        ns = 's' if lat < 0 else 'n'
1✔
1702
        for lon in lons:
1✔
1703
            # east or west?
1704
            ew = 'w' if lon < 0 else 'e'
1✔
1705
            filename = '{}{:02.0f}{}{:03.0f}'.format(ns, abs(lat), ew,
1✔
1706
                                                     abs(lon))
1707
            zones.append(filename)
1✔
1708
    return list(sorted(set(zones)))
1✔
1709

1710

1711
def mapzen_zone(lon_ex, lat_ex, dx_meter=None, zoom=None):
9✔
1712
    """Returns a list of AWS mapzen zones covering the desired extent.
1713

1714
    For mapzen one has to specify the level of detail (zoom) one wants. The
1715
    best way in OGGM is to specify dx_meter of the underlying map and OGGM
1716
    will decide which zoom level works best.
1717
    """
1718

1719
    if dx_meter is None and zoom is None:
1✔
1720
        raise InvalidParamsError('Need either zoom level or dx_meter.')
×
1721

1722
    bottom, top = lat_ex
1✔
1723
    left, right = lon_ex
1✔
1724
    ybound = 85.0511
1✔
1725
    if bottom <= -ybound:
1✔
1726
        bottom = -ybound
1✔
1727
    if top <= -ybound:
1✔
1728
        top = -ybound
1✔
1729
    if bottom > ybound:
1✔
1730
        bottom = ybound
1✔
1731
    if top > ybound:
1✔
1732
        top = ybound
1✔
1733
    if right >= 180:
1✔
1734
        right = 179.999
1✔
1735
    if left >= 180:
1✔
1736
        left = 179.999
1✔
1737

1738
    if dx_meter:
1✔
1739
        # Find out the zoom so that we are close to the desired accuracy
1740
        lat = np.max(np.abs([bottom, top]))
1✔
1741
        zoom = int(np.ceil(math.log2((math.cos(lat * math.pi / 180) *
1✔
1742
                                      2 * math.pi * WEB_EARTH_RADUIS) /
1743
                                     (WEB_N_PIX * dx_meter))))
1744

1745
        # According to this we should just always stay above 10 (sorry)
1746
        # https://github.com/tilezen/joerd/blob/master/docs/data-sources.md
1747
        zoom = 10 if zoom < 10 else zoom
1✔
1748

1749
    # Code from planetutils
1750
    size = 2 ** zoom
1✔
1751
    xt = lambda x: int((x + 180.0) / 360.0 * size)
1✔
1752
    yt = lambda y: int((1.0 - math.log(math.tan(math.radians(y)) +
1✔
1753
                                       (1 / math.cos(math.radians(y))))
1754
                        / math.pi) / 2.0 * size)
1755
    tiles = []
1✔
1756
    for x in range(xt(left), xt(right) + 1):
1✔
1757
        for y in range(yt(top), yt(bottom) + 1):
1✔
1758
            tiles.append('/'.join(map(str, [zoom, x, str(y) + '.tif'])))
1✔
1759
    return tiles
1✔
1760

1761

1762
def get_demo_file(fname):
9✔
1763
    """Returns the path to the desired OGGM-sample-file.
1764

1765
    If Sample data is not cached it will be downloaded from
1766
    https://github.com/OGGM/oggm-sample-data
1767

1768
    Parameters
1769
    ----------
1770
    fname : str
1771
        Filename of the desired OGGM-sample-file
1772

1773
    Returns
1774
    -------
1775
    str
1776
        Absolute path to the desired file.
1777
    """
1778

1779
    d = download_oggm_files()
7✔
1780
    if fname in d:
7✔
1781
        return d[fname]
7✔
1782
    else:
1783
        return None
1✔
1784

1785

1786
def get_wgms_files():
9✔
1787
    """Get the path to the default WGMS-RGI link file and the data dir.
1788

1789
    Returns
1790
    -------
1791
    (file, dir) : paths to the files
1792
    """
1793

1794
    download_oggm_files()
4✔
1795
    sdir = os.path.join(cfg.CACHE_DIR,
4✔
1796
                        'oggm-sample-data-%s' % SAMPLE_DATA_COMMIT,
1797
                        'wgms')
1798
    datadir = os.path.join(sdir, 'mbdata')
4✔
1799
    assert os.path.exists(datadir)
4✔
1800

1801
    outf = os.path.join(sdir, 'rgi_wgms_links_20220112.csv')
4✔
1802
    outf = pd.read_csv(outf, dtype={'RGI_REG': object})
4✔
1803

1804
    return outf, datadir
4✔
1805

1806

1807
def get_glathida_file():
9✔
1808
    """Get the path to the default GlaThiDa-RGI link file.
1809

1810
    Returns
1811
    -------
1812
    file : paths to the file
1813
    """
1814

1815
    # Roll our own
1816
    download_oggm_files()
1✔
1817
    sdir = os.path.join(cfg.CACHE_DIR,
1✔
1818
                        'oggm-sample-data-%s' % SAMPLE_DATA_COMMIT,
1819
                        'glathida')
1820
    outf = os.path.join(sdir, 'rgi_glathida_links.csv')
1✔
1821
    assert os.path.exists(outf)
1✔
1822
    return outf
1✔
1823

1824

1825
def get_rgi_dir(version=None, reset=False):
9✔
1826
    """Path to the RGI directory.
1827

1828
    If the RGI files are not present, download them.
1829

1830
    Parameters
1831
    ----------
1832
    version : str
1833
        '5', '6', defaults to None (linking to the one specified in cfg.PARAMS)
1834
    reset : bool
1835
        If True, deletes the RGI directory first and downloads the data
1836

1837
    Returns
1838
    -------
1839
    str
1840
        path to the RGI directory
1841
    """
1842

1843
    with get_lock():
1✔
1844
        return _get_rgi_dir_unlocked(version=version, reset=reset)
1✔
1845

1846

1847
def _get_rgi_dir_unlocked(version=None, reset=False):
9✔
1848

1849
    rgi_dir = cfg.PATHS['rgi_dir']
1✔
1850
    if version is None:
1✔
1851
        version = cfg.PARAMS['rgi_version']
×
1852

1853
    if len(version) == 1:
1✔
1854
        version += '0'
1✔
1855

1856
    # Be sure the user gave a sensible path to the RGI dir
1857
    if not rgi_dir:
1✔
1858
        raise InvalidParamsError('The RGI data directory has to be'
×
1859
                                 'specified explicitly.')
1860
    rgi_dir = os.path.abspath(os.path.expanduser(rgi_dir))
1✔
1861
    rgi_dir = os.path.join(rgi_dir, 'RGIV' + version)
1✔
1862
    mkdir(rgi_dir, reset=reset)
1✔
1863

1864
    if version == '50':
1✔
1865
        dfile = 'http://www.glims.org/RGI/rgi50_files/rgi50.zip'
1✔
1866
    elif version == '60':
1✔
1867
        dfile = 'http://www.glims.org/RGI/rgi60_files/00_rgi60.zip'
1✔
1868
    elif version == '61':
×
1869
        dfile = 'https://cluster.klima.uni-bremen.de/data/rgi/rgi_61.zip'
×
1870
    elif version == '62':
×
1871
        dfile = 'https://cluster.klima.uni-bremen.de/~oggm/rgi/rgi62.zip'
×
1872

1873
    test_file = os.path.join(rgi_dir,
1✔
1874
                             '*_rgi*{}_manifest.txt'.format(version))
1875

1876
    if len(glob.glob(test_file)) == 0:
1✔
1877
        # if not there download it
1878
        ofile = file_downloader(dfile, reset=reset)
1✔
1879
        # Extract root
1880
        with zipfile.ZipFile(ofile) as zf:
1✔
1881
            zf.extractall(rgi_dir)
1✔
1882
        # Extract subdirs
1883
        pattern = '*_rgi{}_*.zip'.format(version)
1✔
1884
        for root, dirs, files in os.walk(cfg.PATHS['rgi_dir']):
1✔
1885
            for filename in fnmatch.filter(files, pattern):
1✔
1886
                zfile = os.path.join(root, filename)
1✔
1887
                with zipfile.ZipFile(zfile) as zf:
1✔
1888
                    ex_root = zfile.replace('.zip', '')
1✔
1889
                    mkdir(ex_root)
1✔
1890
                    zf.extractall(ex_root)
1✔
1891
                # delete the zipfile after success
1892
                os.remove(zfile)
1✔
1893
        if len(glob.glob(test_file)) == 0:
1✔
1894
            raise RuntimeError('Could not find a manifest file in the RGI '
×
1895
                               'directory: ' + rgi_dir)
1896
    return rgi_dir
1✔
1897

1898

1899
def get_rgi_region_file(region, version=None, reset=False):
9✔
1900
    """Path to the RGI region file.
1901

1902
    If the RGI files are not present, download them.
1903

1904
    Parameters
1905
    ----------
1906
    region : str
1907
        from '01' to '19'
1908
    version : str
1909
        '5', '6', defaults to None (linking to the one specified in cfg.PARAMS)
1910
    reset : bool
1911
        If True, deletes the RGI directory first and downloads the data
1912

1913
    Returns
1914
    -------
1915
    str
1916
        path to the RGI shapefile
1917
    """
1918

1919
    rgi_dir = get_rgi_dir(version=version, reset=reset)
1✔
1920
    f = list(glob.glob(rgi_dir + "/*/*{}_*.shp".format(region)))
1✔
1921
    assert len(f) == 1
1✔
1922
    return f[0]
1✔
1923

1924

1925
def get_rgi_glacier_entities(rgi_ids, version=None):
9✔
1926
    """Get a list of glacier outlines selected from their RGI IDs.
1927

1928
    Will download RGI data if not present.
1929

1930
    Parameters
1931
    ----------
1932
    rgi_ids : list of str
1933
        the glaciers you want the outlines for
1934
    version : str
1935
        the rgi version
1936

1937
    Returns
1938
    -------
1939
    geopandas.GeoDataFrame
1940
        containing the desired RGI glacier outlines
1941
    """
1942

1943
    regions = [s.split('-')[1].split('.')[0] for s in rgi_ids]
×
1944
    if version is None:
×
1945
        version = rgi_ids[0].split('-')[0][-2:]
×
1946
    selection = []
×
1947
    for reg in sorted(np.unique(regions)):
×
1948
        sh = gpd.read_file(get_rgi_region_file(reg, version=version))
×
1949
        selection.append(sh.loc[sh.RGIId.isin(rgi_ids)])
×
1950

1951
    # Make a new dataframe of those
1952
    selection = pd.concat(selection)
×
1953
    selection.crs = sh.crs  # for geolocalisation
×
1954
    if len(selection) != len(rgi_ids):
×
1955
        raise RuntimeError('Could not find all RGI ids')
×
1956

1957
    return selection
×
1958

1959

1960
def get_rgi_intersects_dir(version=None, reset=False):
9✔
1961
    """Path to the RGI directory containing the intersect files.
1962

1963
    If the files are not present, download them.
1964

1965
    Parameters
1966
    ----------
1967
    version : str
1968
        '5', '6', defaults to None (linking to the one specified in cfg.PARAMS)
1969
    reset : bool
1970
        If True, deletes the intersects before redownloading them
1971

1972
    Returns
1973
    -------
1974
    str
1975
        path to the directory
1976
    """
1977

1978
    with get_lock():
1✔
1979
        return _get_rgi_intersects_dir_unlocked(version=version, reset=reset)
1✔
1980

1981

1982
def _get_rgi_intersects_dir_unlocked(version=None, reset=False):
9✔
1983

1984
    rgi_dir = cfg.PATHS['rgi_dir']
1✔
1985
    if version is None:
1✔
1986
        version = cfg.PARAMS['rgi_version']
×
1987

1988
    if len(version) == 1:
1✔
1989
        version += '0'
1✔
1990

1991
    # Be sure the user gave a sensible path to the RGI dir
1992
    if not rgi_dir:
1✔
1993
        raise InvalidParamsError('The RGI data directory has to be'
×
1994
                                 'specified explicitly.')
1995

1996
    rgi_dir = os.path.abspath(os.path.expanduser(rgi_dir))
1✔
1997
    mkdir(rgi_dir)
1✔
1998

1999
    dfile = 'https://cluster.klima.uni-bremen.de/data/rgi/'
1✔
2000
    dfile += 'RGI_V{}_Intersects.zip'.format(version)
1✔
2001
    if version == '62':
1✔
2002
        dfile = ('https://cluster.klima.uni-bremen.de/~oggm/rgi/'
×
2003
                 'rgi62_Intersects.zip')
2004

2005
    odir = os.path.join(rgi_dir, 'RGI_V' + version + '_Intersects')
1✔
2006
    if reset and os.path.exists(odir):
1✔
2007
        shutil.rmtree(odir)
×
2008

2009
    # A lot of code for backwards compat (sigh...)
2010
    if version in ['50', '60']:
1✔
2011
        test_file = os.path.join(odir, 'Intersects_OGGM_Manifest.txt')
1✔
2012
        if not os.path.exists(test_file):
1✔
2013
            # if not there download it
2014
            ofile = file_downloader(dfile, reset=reset)
1✔
2015
            # Extract root
2016
            with zipfile.ZipFile(ofile) as zf:
1✔
2017
                zf.extractall(odir)
1✔
2018
            if not os.path.exists(test_file):
1✔
2019
                raise RuntimeError('Could not find a manifest file in the RGI '
×
2020
                                   'directory: ' + odir)
2021
    else:
2022
        test_file = os.path.join(odir,
×
2023
                                 '*ntersect*anifest.txt'.format(version))
2024
        if len(glob.glob(test_file)) == 0:
×
2025
            # if not there download it
2026
            ofile = file_downloader(dfile, reset=reset)
×
2027
            # Extract root
2028
            with zipfile.ZipFile(ofile) as zf:
×
2029
                zf.extractall(odir)
×
2030
            # Extract subdirs
2031
            pattern = '*_rgi{}_*.zip'.format(version)
×
2032
            for root, dirs, files in os.walk(cfg.PATHS['rgi_dir']):
×
2033
                for filename in fnmatch.filter(files, pattern):
×
2034
                    zfile = os.path.join(root, filename)
×
2035
                    with zipfile.ZipFile(zfile) as zf:
×
2036
                        ex_root = zfile.replace('.zip', '')
×
2037
                        mkdir(ex_root)
×
2038
                        zf.extractall(ex_root)
×
2039
                    # delete the zipfile after success
2040
                    os.remove(zfile)
×
2041
            if len(glob.glob(test_file)) == 0:
×
2042
                raise RuntimeError('Could not find a manifest file in the RGI '
×
2043
                                   'directory: ' + odir)
2044

2045
    return odir
1✔
2046

2047

2048
def get_rgi_intersects_region_file(region=None, version=None, reset=False):
9✔
2049
    """Path to the RGI regional intersect file.
2050

2051
    If the RGI files are not present, download them.
2052

2053
    Parameters
2054
    ----------
2055
    region : str
2056
        from '00' to '19', with '00' being the global file (deprecated).
2057
        From RGI version '61' onwards, please use `get_rgi_intersects_entities`
2058
        with a list of glaciers instead of relying to the global file.
2059
    version : str
2060
        '5', '6', '61'... defaults the one specified in cfg.PARAMS
2061
    reset : bool
2062
        If True, deletes the intersect file before redownloading it
2063

2064
    Returns
2065
    -------
2066
    str
2067
        path to the RGI intersects shapefile
2068
    """
2069

2070
    if version is None:
1✔
2071
        version = cfg.PARAMS['rgi_version']
×
2072
    if len(version) == 1:
1✔
2073
        version += '0'
1✔
2074

2075
    rgi_dir = get_rgi_intersects_dir(version=version, reset=reset)
1✔
2076

2077
    if region == '00':
1✔
2078
        if version in ['50', '60']:
1✔
2079
            version = 'AllRegs'
1✔
2080
            region = '*'
1✔
2081
        else:
2082
            raise InvalidParamsError("From RGI version 61 onwards, please use "
×
2083
                                     "get_rgi_intersects_entities() instead.")
2084
    f = list(glob.glob(os.path.join(rgi_dir, "*", '*intersects*' + region +
1✔
2085
                                    '_rgi*' + version + '*.shp')))
2086
    assert len(f) == 1
1✔
2087
    return f[0]
1✔
2088

2089

2090
def get_rgi_intersects_entities(rgi_ids, version=None):
9✔
2091
    """Get a list of glacier intersects selected from their RGI IDs.
2092

2093
    Parameters
2094
    ----------
2095
    rgi_ids: list of str
2096
        list of rgi_ids you want to look for intersections for
2097
    version: str
2098
        '5', '6', '61'... defaults the one specified in cfg.PARAMS
2099

2100
    Returns
2101
    -------
2102
    geopandas.GeoDataFrame
2103
        with the selected intersects
2104
    """
2105

2106
    if version is None:
×
2107
        version = cfg.PARAMS['rgi_version']
×
2108
    if len(version) == 1:
×
2109
        version += '0'
×
2110
    try:
×
2111
        regions = [s.split('-')[3] for s in rgi_ids]
×
2112

2113
    except IndexError:
×
2114
        # RGI V6
2115
        regions = [s.split('-')[1].split('.')[0] for s in rgi_ids]
×
2116
    selection = []
×
2117
    for reg in sorted(np.unique(regions)):
×
2118
        sh = gpd.read_file(get_rgi_intersects_region_file(reg,
×
2119
                                                          version=version))
2120
        selection.append(sh.loc[sh.RGIId_1.isin(rgi_ids) |
×
2121
                                sh.RGIId_2.isin(rgi_ids)])
2122

2123
    # Make a new dataframe of those
2124
    selection = pd.concat(selection)
×
2125
    selection.crs = sh.crs  # for geolocalisation
×
2126

2127
    return selection
×
2128

2129

2130
def is_dem_source_available(source, lon_ex, lat_ex):
9✔
2131
    """Checks if a DEM source is available for your purpose.
2132

2133
    This is only a very rough check! It doesn't mean that the data really is
2134
    available, but at least it's worth a try.
2135

2136
    Parameters
2137
    ----------
2138
    source : str, required
2139
        the source you want to check for
2140
    lon_ex : tuple or int, required
2141
        a (min_lon, max_lon) tuple delimiting the requested area longitudes
2142
    lat_ex : tuple or int, required
2143
        a (min_lat, max_lat) tuple delimiting the requested area latitudes
2144

2145
    Returns
2146
    -------
2147
    True or False
2148
    """
2149
    from oggm.utils import tolist
6✔
2150
    lon_ex = tolist(lon_ex, length=2)
6✔
2151
    lat_ex = tolist(lat_ex, length=2)
6✔
2152

2153
    def _in_grid(grid_json, lon, lat):
6✔
2154
        i, j = cfg.DATA['dem_grids'][grid_json].transform(lon, lat,
1✔
2155
                                                          maskout=True)
2156
        return np.all(~ (i.mask | j.mask))
1✔
2157

2158
    if source == 'GIMP':
6✔
2159
        return _in_grid('gimpdem_90m_v01.1.json', lon_ex, lat_ex)
1✔
2160
    elif source == 'ARCTICDEM':
6✔
2161
        return _in_grid('arcticdem_mosaic_100m_v3.0.json', lon_ex, lat_ex)
1✔
2162
    elif source == 'RAMP':
6✔
2163
        return _in_grid('AntarcticDEM_wgs84.json', lon_ex, lat_ex)
1✔
2164
    elif source == 'REMA':
6✔
2165
        return _in_grid('REMA_100m_dem.json', lon_ex, lat_ex)
1✔
2166
    elif source == 'ALASKA':
6✔
2167
        return _in_grid('Alaska_albers_V3.json', lon_ex, lat_ex)
×
2168
    elif source == 'TANDEM':
6✔
2169
        return True
1✔
2170
    elif source == 'AW3D30':
6✔
2171
        return np.min(lat_ex) > -60
1✔
2172
    elif source == 'MAPZEN':
6✔
2173
        return True
1✔
2174
    elif source == 'DEM3':
6✔
2175
        return True
1✔
2176
    elif source == 'ASTER':
6✔
2177
        return True
1✔
2178
    elif source == 'SRTM':
6✔
2179
        return np.max(np.abs(lat_ex)) < 60
1✔
2180
    elif source in ['COPDEM30', 'COPDEM90']:
6✔
2181
        return True
×
2182
    elif source == 'NASADEM':
6✔
2183
        return (np.min(lat_ex) > -56) and (np.max(lat_ex) < 60)
×
2184
    elif source == 'USER':
6✔
UNCOV
2185
        return True
×
2186
    elif source is None:
6✔
2187
        return True
6✔
2188

2189

2190
def default_dem_source(rgi_id):
9✔
2191
    """Current default DEM source at a given location.
2192

2193
    Parameters
2194
    ----------
2195
    rgi_id : str
2196
        the RGI id
2197

2198
    Returns
2199
    -------
2200
    the chosen DEM source
2201
    """
2202
    rgi_reg = 'RGI{}'.format(rgi_id[6:8])
1✔
2203
    rgi_id = rgi_id[:14]
1✔
2204
    if cfg.DEM_SOURCE_TABLE.get(rgi_reg) is None:
1✔
2205
        fp = get_demo_file('rgi62_dem_frac.h5')
1✔
2206
        cfg.DEM_SOURCE_TABLE[rgi_reg] = pd.read_hdf(fp)
1✔
2207

2208
    sel = cfg.DEM_SOURCE_TABLE[rgi_reg].loc[rgi_id]
1✔
2209
    for s in ['NASADEM', 'COPDEM90', 'COPDEM30', 'GIMP', 'REMA',
1✔
2210
              'RAMP', 'TANDEM', 'MAPZEN']:
2211
        if sel.loc[s] > 0.75:
1✔
2212
            return s
1✔
2213
    # If nothing works, try COPDEM again
2214
    return 'COPDEM90'
×
2215

2216

2217
def get_topo_file(lon_ex=None, lat_ex=None, rgi_id=None, *,
9✔
2218
                  dx_meter=None, zoom=None, source=None):
2219
    """Path(s) to the DEM file(s) covering the desired extent.
2220

2221
    If the needed files for covering the extent are not present, download them.
2222

2223
    The default behavior is to try a list of DEM sources in order, and
2224
    stop once the downloaded data is covering a large enough part of the
2225
    glacier. The DEM sources are tested in the following order:
2226

2227
    'NASADEM' -> 'COPDEM' -> 'GIMP' -> 'REMA' -> 'TANDEM' -> 'MAPZEN'
2228

2229
    To force usage of a certain data source, use the ``source`` kwarg argument.
2230

2231
    Parameters
2232
    ----------
2233
    lon_ex : tuple or int, required
2234
        a (min_lon, max_lon) tuple delimiting the requested area longitudes
2235
    lat_ex : tuple or int, required
2236
        a (min_lat, max_lat) tuple delimiting the requested area latitudes
2237
    rgi_id : str, required if source=None
2238
        the glacier id, used to decide on the DEM source
2239
    rgi_region : str, optional
2240
        the RGI region number (required for the GIMP DEM)
2241
    rgi_subregion : str, optional
2242
        the RGI subregion str (useful for RGI Reg 19)
2243
    dx_meter : float, required for source='MAPZEN'
2244
        the resolution of the glacier map (to decide the zoom level of mapzen)
2245
    zoom : int, optional
2246
        if you know the zoom already (for MAPZEN only)
2247
    source : str or list of str, optional
2248
        Name of specific DEM source. see utils.DEM_SOURCES for a list
2249

2250
    Returns
2251
    -------
2252
    tuple: (list with path(s) to the DEM file(s), data source str)
2253
    """
2254
    from oggm.utils import tolist
6✔
2255
    lon_ex = tolist(lon_ex, length=2)
6✔
2256
    lat_ex = tolist(lat_ex, length=2)
6✔
2257

2258
    if source is not None and not isinstance(source, str):
6✔
2259
        # check all user options
2260
        for s in source:
×
2261
            demf, source_str = get_topo_file(lon_ex=lon_ex, lat_ex=lat_ex,
×
2262
                                             rgi_id=rgi_id, source=s)
2263
            if demf[0]:
×
2264
                return demf, source_str
×
2265

2266
    # Did the user specify a specific DEM file?
2267
    if 'dem_file' in cfg.PATHS and os.path.isfile(cfg.PATHS['dem_file']):
6✔
2268
        source = 'USER' if source is None else source
6✔
2269
        if source == 'USER':
6✔
2270
            return [cfg.PATHS['dem_file']], source
6✔
2271

2272
    # Some logic to decide which source to take if unspecified
2273
    if source is None:
1✔
2274
        if rgi_id is None:
×
2275
            raise InvalidParamsError('rgi_id is needed if source=None')
×
2276
        source = default_dem_source(rgi_id)
×
2277

2278
    if source not in DEM_SOURCES:
1✔
2279
        raise InvalidParamsError('`source` must be one of '
×
2280
                                 '{}'.format(DEM_SOURCES))
2281

2282
    # OK go
2283
    files = []
1✔
2284
    if source == 'GIMP':
1✔
2285
        _file = _download_topo_file_from_cluster('gimpdem_90m_v01.1.tif')
1✔
2286
        files.append(_file)
1✔
2287

2288
    if source == 'ARCTICDEM':
1✔
2289
        zones = arcticdem_zone(lon_ex, lat_ex)
1✔
2290
        for z in zones:
1✔
2291
            with get_lock():
1✔
2292
                url = 'https://cluster.klima.uni-bremen.de/~oggm/'
1✔
2293
                url += 'dem/ArcticDEM_100m_v3.0/'
1✔
2294
                url += '{}_100m_v3.0/{}_100m_v3.0_reg_dem.tif'.format(z, z)
1✔
2295
                files.append(file_downloader(url))
1✔
2296

2297
    if source == 'RAMP':
1✔
2298
        _file = _download_topo_file_from_cluster('AntarcticDEM_wgs84.tif')
1✔
2299
        files.append(_file)
1✔
2300

2301
    if source == 'ALASKA':
1✔
2302
        zones = alaska_dem_zone(lon_ex, lat_ex)
1✔
2303
        for z in zones:
1✔
2304
            with get_lock():
1✔
2305
                url = 'https://cluster.klima.uni-bremen.de/~oggm/'
1✔
2306
                url += 'dem/Alaska_albers_V3/'
1✔
2307
                url += '{}_Alaska_albers_V3/'.format(z)
1✔
2308
                url += '{}_Alaska_albers_V3.tif'.format(z)
1✔
2309
                files.append(file_downloader(url))
1✔
2310

2311
    if source == 'REMA':
1✔
2312
        zones = rema_zone(lon_ex, lat_ex)
1✔
2313
        for z in zones:
1✔
2314
            with get_lock():
1✔
2315
                url = 'https://cluster.klima.uni-bremen.de/~oggm/'
1✔
2316
                url += 'dem/REMA_100m_v1.1/'
1✔
2317
                url += '{}_100m_v1.1/{}_100m_v1.1_reg_dem.tif'.format(z, z)
1✔
2318
                files.append(file_downloader(url))
1✔
2319

2320
    if source == 'TANDEM':
1✔
2321
        zones = tandem_zone(lon_ex, lat_ex)
1✔
2322
        for z in zones:
1✔
2323
            files.append(_download_tandem_file(z))
1✔
2324

2325
    if source == 'AW3D30':
1✔
2326
        zones = aw3d30_zone(lon_ex, lat_ex)
1✔
2327
        for z in zones:
1✔
2328
            files.append(_download_aw3d30_file(z))
1✔
2329

2330
    if source == 'MAPZEN':
1✔
2331
        zones = mapzen_zone(lon_ex, lat_ex, dx_meter=dx_meter, zoom=zoom)
1✔
2332
        for z in zones:
1✔
2333
            files.append(_download_mapzen_file(z))
1✔
2334

2335
    if source == 'ASTER':
1✔
2336
        zones = aster_zone(lon_ex, lat_ex)
1✔
2337
        for z in zones:
1✔
2338
            files.append(_download_aster_file(z))
1✔
2339

2340
    if source == 'DEM3':
1✔
2341
        zones = dem3_viewpano_zone(lon_ex, lat_ex)
1✔
2342
        for z in zones:
1✔
2343
            files.append(_download_dem3_viewpano(z))
1✔
2344

2345
    if source == 'SRTM':
1✔
2346
        zones = srtm_zone(lon_ex, lat_ex)
1✔
2347
        for z in zones:
1✔
2348
            files.append(_download_srtm_file(z))
1✔
2349

2350
    if source in ['COPDEM30', 'COPDEM90']:
1✔
2351
        filetuple = copdem_zone(lon_ex, lat_ex, source)
1✔
2352
        for cpp, eop in filetuple:
1✔
2353
            files.append(_download_copdem_file(cpp, eop, source))
1✔
2354

2355
    if source == 'NASADEM':
1✔
2356
        zones = nasadem_zone(lon_ex, lat_ex)
1✔
2357
        for z in zones:
1✔
2358
            files.append(_download_nasadem_file(z))
1✔
2359

2360
    # filter for None (e.g. oceans)
2361
    files = [s for s in files if s]
1✔
2362
    if files:
1✔
2363
        return files, source
1✔
2364
    else:
2365
        raise InvalidDEMError('Source: {2} no topography file available for '
×
2366
                              'extent lat:{0}, lon:{1}!'.
2367
                              format(lat_ex, lon_ex, source))
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc