Source code for oggm.shop.histalp

import logging
import warnings

# External libs
import numpy as np
import pandas as pd
from scipy import stats

# Optional libs
try:
    import salem
except ImportError:
    pass

# Locals
from oggm import cfg
from oggm import utils
from oggm import entity_task
from oggm.exceptions import InvalidParamsError

# Module logger
log = logging.getLogger(__name__)

HISTALP_SERVER = 'http://www.zamg.ac.at/histalp/download/grid5m/'


def set_histalp_url(url):
    """If you want to use a different server for HISTALP (for testing, etc)."""
    global HISTALP_SERVER
    HISTALP_SERVER = url


[docs]@utils.locked_func def get_histalp_file(var=None): """Returns a path to the desired HISTALP baseline climate file. If the file is not present, download it. Parameters ---------- var : str 'tmp' for temperature 'pre' for precipitation Returns ------- str path to the file """ # Be sure input makes sense if var not in ['tmp', 'pre']: raise InvalidParamsError('HISTALP variable {} ' 'does not exist!'.format(var)) # File to look for if var == 'tmp': bname = 'HISTALP_temperature_1780-2014.nc' else: bname = 'HISTALP_precipitation_all_abs_1801-2014.nc' h_url = HISTALP_SERVER + bname + '.bz2' return utils.file_extractor(utils.file_downloader(h_url))
[docs]@entity_task(log, writes=['climate_historical']) def process_histalp_data(gdir, y0=None, y1=None, output_filesuffix=None): """Processes and writes the HISTALP baseline climate data for this glacier. Extracts the nearest timeseries and writes everything to a NetCDF file. Parameters ---------- gdir : :py:class:`oggm.GlacierDirectory` the glacier directory to process y0 : int the starting year of the timeseries to write. The default is to take 1850 (because the data is quite bad before that) y1 : int the ending year of the timeseries to write. The default is to take the entire time period available in the file, but with this kwarg you can shorten it (to save space or to crop bad data) output_filesuffix : str this adds a suffix to the output file (useful to avoid overwriting previous experiments) """ if cfg.PARAMS['baseline_climate'] != 'HISTALP': raise InvalidParamsError("cfg.PARAMS['baseline_climate'] should be " "set to HISTALP.") # read the time out of the pure netcdf file ft = get_histalp_file('tmp') fp = get_histalp_file('pre') with utils.ncDataset(ft) as nc: vt = nc.variables['time'] assert vt[0] == 0 assert vt[-1] == vt.shape[0] - 1 t0 = vt.units.split(' since ')[1][:7] time_t = pd.date_range(start=t0, periods=vt.shape[0], freq='MS') with utils.ncDataset(fp) as nc: vt = nc.variables['time'] assert vt[0] == 0.5 assert vt[-1] == vt.shape[0] - .5 t0 = vt.units.split(' since ')[1][:7] time_p = pd.date_range(start=t0, periods=vt.shape[0], freq='MS') # Now open with salem nc_ts_tmp = salem.GeoNetcdf(ft, time=time_t) nc_ts_pre = salem.GeoNetcdf(fp, time=time_p) # Some default if y0 is None: y0 = 1850 # set temporal subset for the ts data (hydro years) # the reference time is given by precip, which is shorter yrs = nc_ts_pre.time.year y0 = yrs[0] if y0 is None else y0 y1 = yrs[-1] if y1 is None else y1 nc_ts_tmp.set_period(t0=f'{y0}-01-01', t1=f'{y1}-12-01') nc_ts_pre.set_period(t0=f'{y0}-01-01', t1=f'{y1}-12-01') time = nc_ts_pre.time ny, r = divmod(len(time), 12) assert r == 0 # Units assert nc_ts_tmp._nc.variables['HSURF'].units.lower() in ['m', 'meters', 'meter', 'metres', 'metre'] assert nc_ts_tmp._nc.variables['T_2M'].units.lower() in ['degc', 'degrees', 'degrees celcius', 'degree', 'c'] assert nc_ts_pre._nc.variables['TOT_PREC'].units.lower() in ['kg m-2', 'l m-2', 'mm', 'millimeters', 'millimeter'] # geoloc lon = gdir.cenlon lat = gdir.cenlat nc_ts_tmp.set_subset(corners=((lon, lat), (lon, lat)), margin=1) nc_ts_pre.set_subset(corners=((lon, lat), (lon, lat)), margin=1) # read the data temp = nc_ts_tmp.get_vardata('T_2M') prcp = nc_ts_pre.get_vardata('TOT_PREC') hgt = nc_ts_tmp.get_vardata('HSURF') ref_lon = nc_ts_tmp.get_vardata('lon') ref_lat = nc_ts_tmp.get_vardata('lat') source = nc_ts_tmp._nc.title[:7] nc_ts_tmp._nc.close() nc_ts_pre._nc.close() gdir.write_monthly_climate_file(time, prcp[:, 1, 1], temp[:, 1, 1], hgt[1, 1], ref_lon[1], ref_lat[1], filesuffix=output_filesuffix, source=source)