Add pyaps3 files

__init__.py

@@ -0,0 +1,21 @@
+'''
+PyAPS module to compute InSAR phase delay maps from weather models.
+
+Written by Romain Jolivet <rjolivet@gps.caltech.edu> and Piyush Agram <piyush@gps.caltech.edu>. The original Fortran package was written by Romain Jolivet and the Python version including support for different models was written by Piyush Agram.
+
+.. note::
+	Details of the python module can be obtained `here. <http://code.google.com/p/pyaps>`_
+'''
+
+__all__ = ['geocoord','rdrcoord','autoget']
+
+from .geocoord import PyAPS_geo
+from .rdrcoord import PyAPS_rdr
+from .autoget import *
+
+
+############################################################
+# Program is part of PyAPS                                 #
+# Copyright 2012, by the California Institute of Technology#
+# Contact: earthdef@gps.caltech.edu                        #
+############################################################

autoget.py

@@ -0,0 +1,149 @@
+############################################################
+# Program is part of PyAPS                                 #
+# Copyright 2012, by the California Institute of Technology#
+# Contact: earthdef@gps.caltech.edu                        #
+# Modified by A. Benoit and R. Jolivet 2019                #
+# Ecole Normale Superieure, Paris                          #
+# Contact: insar@geologie.ens.fr                           #
+############################################################
+
+import configparser as ConfigParser
+import sys
+import os.path
+import cdsapi
+
+def ECMWFdload(bdate,hr,filedir,model,datatype,humidity='Q'):
+    '''
+    ECMWF data downloading.
+
+    Args:
+        * bdate     : date to download (str)
+        * hr        : hour to download (str)
+        * filedir   : files directory (str)
+        * model     : weather model ('interim' or 'era5' or 'hres')
+        * datatype  : reanalysis data type (an)
+        * humidity  : humidity
+    '''
+
+    #-------------------------------------------
+    # Initialize
+
+    # Check data
+    assert model in ('era5','interim','hres'), 'Unknown model for ECMWF'
+
+    # Infos for downloading
+    if model in 'interim':
+        print('WARNING: you are downloading from the old ECMWF platform. ERA-Interim is deprecated, use ERA-5 instead.')
+    if model in 'era5':
+        print('INFO: You are using the latest ECMWF platform for downloading datasets: https://cds.climate.copernicus.eu/api/v2')
+
+    #-------------------------------------------
+    # Define parameters
+
+    # Humidity
+    assert humidity in ('Q','R'), 'Unknown humidity field for ECMWF'
+    if humidity in 'Q':
+        if model in 'era5':
+            humidparam = 'specific_humidity'
+        else:
+            humidparam = 133
+    elif humidity in 'R':
+        if model in 'era5':
+            humidparam = 'relative_humidity'
+        else:
+            humidparam = 157
+
+    # Grid size (only for HRES and ERA-Interim)
+    if model in 'hres':
+        gridsize = '0.10/0.10'
+    elif model in 'interim':
+        gridsize = '0.75/0.75'
+
+    #-------------------------------------------
+    # Iterate over dates
+
+    flist = []
+    for k in range(len(bdate)):
+        day = bdate[k]
+
+        #-------------------------------------------
+        # CASE 1: request for CDS API client (new ECMWF platform, for ERA-5)    
+        if model in 'era5':
+
+            # Contact the server
+            c = cdsapi.Client()
+
+            # Pressure levels
+            pressure_lvls = ['1','2','3','5','7','10','20','30','50', 
+                            '70','100','125','150','175','200','225',
+                            '250','300','350','400','450','500','550',
+                            '600','650','700','750','775','800','825',
+                            '850','875','900','925','950','975','1000']
+
+            # Dictionary
+            indict = {'product_type':'reanalysis',
+                      'format':'grib',
+                      'variable':['geopotential','temperature','{}'.format(humidparam)],
+                      'pressure_level': pressure_lvls,
+                      'year':'{}'.format(day[0:4]),
+                      'month':'{}'.format(day[4:6]),
+                      'day':'{}'.format(day[6:8]),
+                      'time':'{}:00'.format(hr)}
+
+            # Output
+            fname = '%s/ERA-5_%s_%s_%s.grb'%(filedir,day,hr,datatype)
+            flist.append(fname)
+
+            # Assert grib file not yet downloaded
+            if not os.path.exists(fname):
+                print('Downloading %d of %d: %s '%(k+1,len(bdate),fname))
+                print(indict)
+
+                # Make the request
+                c.retrieve(
+                        'reanalysis-{}-pressure-levels'.format(model),
+                        indict,
+                        fname)
+
+        #-------------------------------------------
+        # CASE 2: request for WEB API client (old ECMWF platform, deprecated, for ERA-Int and HRES)
+        else:
+
+            # Contact the server
+            from pyaps3.ecmwfapi import ECMWFDataServer
+            dpath = os.path.dirname(__file__)
+            config = ConfigParser.RawConfigParser()
+            config.read(os.path.expanduser('~/.ecmwfcfg'))
+            url = "https://api.ecmwf.int/v1"
+            emid = config.get('ECMWF', 'email')
+            key = config.get('ECMWF', 'key')
+            server = ECMWFDataServer(url=url, key=key, email=emid)
+
+            # Output
+            if model in 'interim':
+                fname = '%s/ERA-Int_%s_%s_%s.grb'%(filedir,day,hr,datatype)
+            elif model in 'hres':
+                fname = '%s/HRES_%s_%s_%s.grb'%(filedir,day,hr,datatype)
+            flist.append(fname)
+
+            # Dictionary
+            indict = {'dataset'  : '{}'.format(model),
+                      'type'     : '{}'.format(datatype),
+                      'date'     : '{}'.format(day),
+                      'time'     : '{}'.format(hr),
+                      'step'     : '0',
+                      'levtype'  : "pl",
+                      'levelist' : "all",
+                      'grid'     : '{}'.format(gridsize),
+                      'param'    : '129/130/{}'.format(humidparam),
+                      'target'   : '{}'.format(fname)}
+
+            # Assert grib file not yet downloaded
+            if not os.path.exists(fname):
+                print('Downloading %d of %d: %s '%(k+1,len(bdate),fname))
+                print(indict)
+
+                # Make the request
+                server.retrieve(indict)
+
+    return flist

delay.py

@@ -0,0 +1,115 @@
+############################################################
+# Program is part of PyAPS                                 #
+# Copyright 2012, by the California Institute of Technology#
+# Contact: earthdef@gps.caltech.edu                        #
+# Modified by A. Benoit and R. Jolivet 2019                #
+# Ecole Normale Superieure, Paris                          #
+# Contact: insar@geologie.ens.fr                           #
+############################################################
+import numpy as np
+
+#############Clausius-Clapeyron for ECMWF as used in Jolivet et al 2011#
+def cc_eraorig(tmp,cdic):
+    '''This routine takes temperature profiles and returns Saturation water vapor partial pressure using the Clausius-Clapeyron law applied in Jolivet et al. 2011,GRL, doi:10.1029/2011GL048757. It can be used in case you are using Relative Humidity from ECMWF models
+
+    Args:
+    	* tmp  (np.array) : Temperature
+    	* cdic (dict)     : Dictionnary of constants
+
+    Returns:
+    	* esat (np.array) : Saturation water vapor partial pressure.'''
+
+    a1w=cdic['a1w']
+    T3=cdic['T3']
+    Rv=cdic['Rv']
+
+    esat=a1w*np.exp( (2.5e6/Rv) * ( (1/T3) - (1/tmp) ) )
+
+    return esat
+###############Completed CC_ERAORIG#####################################
+
+##############Read Input text file #####################################
+def read_eratxt(fname,cdic):
+    '''Read ECMWF files from 0.75 degree grid similar to Romain Jolivet's Delay Package.
+
+    Args:
+        * fname       (str):  Path to the delay file
+        * cdic   (np.float):  Dictionary of constants
+
+    Kwargs:
+        * humidity    (str): Specific ('Q') or relative humidity ('R').
+
+    Returns:
+        * lvls   (np.array): Pressure levels
+        * latlist(np.array): Latitudes of the stations
+        * lonlist(np.array): Longitudes of the stations
+        * gph    (np.array): Geopotential height
+        * tmp    (np.array): Temperature
+        * vpr    (np.array): Vapor pressure
+
+    .. note::
+        Uses cc_eraorig by default.
+    '''
+
+    lvls=[]
+    latlist=[]
+    lonlist=[]
+    gpht=[]
+    tmpt=[]
+    reht=[]
+
+    g=cdic['g']
+
+    f=open(fname,'r')
+    tmp=f.readlines()
+    i=0
+    nstn=0
+    maxloop=np.int(len(tmp))
+    while i<maxloop:
+        if (tmp[i][0]=='-'):
+            nstn=nstn+1
+            lonlat=tmp[i+3].rsplit(' ')
+            lonlist.append(np.float(lonlat[3]))
+            latlist.append(np.float(lonlat[9]))
+            i=i+5
+            new='y'
+        else:
+            if new in ('y'):
+                n=1
+                val=tmp[i].rsplit(' ')
+                gpht.append(np.float(val[0]))
+                lvls.append(np.float(val[1]))
+                tmpt.append(np.float(val[2]))
+                reht.append(np.float(val[3]))
+                i=i+1
+                new='n'
+            else:
+                n=n+1
+                val=tmp[i].rsplit(' ')
+                gpht.append(np.float(val[0]))
+                lvls.append(np.float(val[1]))
+                tmpt.append(np.float(val[2]))
+                reht.append(np.float(val[3]))
+                i=i+1
+
+    gpht=np.array(gpht)/g
+    gph=np.flipud(gpht.reshape((n,nstn),order='F'))
+    del gpht
+
+    tmpt=np.array(tmpt)
+    esat=cc_eraorig(tmpt,cdic)
+    tmp=np.flipud(tmpt.reshape((n,nstn),order='F'))
+    del tmpt
+
+    vprt=(np.array(reht)/100.)*esat
+    vpr=np.flipud(vprt.reshape((n,nstn),order='F'))
+    del vprt
+    del esat    
+
+    lvls=np.flipud(np.array(lvls))
+    lvls=lvls[0:n]
+
+    lonlist=np.array(lonlist)
+    latlist=np.array(latlist)
+
+    return lvls,latlist,lonlist,gph,tmp,vpr

ecmwf.py

@@ -0,0 +1,109 @@
+############################################################
+# Program is part of PyAPS                                 #
+# Copyright 2012, by the California Institute of Technology#
+# Contact: earthdef@gps.caltech.edu                        #
+# Modified by A. Benoit and R. Jolivet 2019                #
+# Ecole Normale Superieure, Paris                          #
+# Contact: insar@geologie.ens.fr                           #
+############################################################
+import urllib.request
+import urllib.parse
+import time
+import datetime
+
+class ECMWFDataServer:
+    def __init__(self,portal,token,email):
+        self.version = '0.3'
+        self.portal  = portal
+        self.token   = token
+        self.email   = email
+
+    def _call(self,action,args):
+
+        params = {'_token'   : self.token,
+                  '_email'   : self.email,
+                  '_action'  : action,
+                  '_version' : self.version}
+        params.update(args)
+
+        data = urllib.parse.urlencode(params)
+        req = urllib.request.Request(self.portal, data)
+        response = urllib.request.urlopen(req)
+
+        json = response.read();
+
+        undef = None;
+        json = eval(json)
+
+        if json != None:
+            if 'error' in json:
+                raise RuntimeError(json['error'])
+            if 'message' in json:
+                self.put(json['message'])
+
+        return json
+
+    def put(self,*args):
+        print(datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')),
+        for a in args:
+            print(a),
+        print()
+
+
+    def retrieve(self,args):
+        self.put("ECMWF data server batch tool version",self.version);
+        user = self._call("user_info",{});
+        self.put("Welcome to",user['name'], "from" , user['organisation']);
+
+        r = self._call('retrieve',args)
+        rid = r['request']
+
+        last  = ''
+        sleep = 0
+        while r['status'] != 'complete' and r['status'] != 'aborted':
+            text = r['status'] + '.'
+            if 'info' in r and r['info'] != None:
+                text = text + ' ' + r['info'] 
+
+            if text != last:
+                self.put("Request",text)
+                last = text
+
+            time.sleep(sleep)
+            r = self._call('status',{'request':rid})
+            if sleep < 60:
+                sleep = sleep + 1
+
+        if r['status'] != last:
+            self.put("Request",r['status'])
+
+        if 'reason' in r:
+            for m in r['reason']:
+                self.put(m)
+
+        if 'result' in r:
+            size = long(r['size'])
+            self.put("Downloading",self._bytename(size))
+            done = self._transfer(r['result'],args['target'])
+            self.put("Done")
+            if done != size:
+                raise RuntimeError("Size mismatch: " + str(done) + " and " + str(size))
+
+        self._call('delete',{'request':rid})
+
+        if r['status'] == 'aborted':
+            raise RuntimeError("Request aborted")
+
+    def _transfer(self,url,path):
+        result =  urllib.request.urlretrieve(url,path)
+        return long(result[1]['content-length'])
+
+    def _bytename(self,size):   
+        next = {'':'K','K':'M','M':'G','G':'T','T':'P'}
+        l    = ''
+        size = size*1.0
+        while 1024 < size:
+            l = next[l]
+            size = size / 1024
+        return "%g %sbyte%s" % (size,l,'s')
+