############################################################ # 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 os.path import numpy as np import scipy.integrate as intg import scipy.interpolate as si import scipy.spatial as ss import pyaps3.processor as processor import pyaps3.utils as utils def sanity_check(model): if model in ('ECMWF','ERA','NARR'): from . import era, narr if model=='MERRA': from . import merra return ##############Creating a class object for PyAPS use. class PyAPS_geo: ''' Class for dealing with Atmospheric phase corrections in geo-coded space. Operates on one weather model file and one Geo.rsc file at a time. This use the relative humidity information from the specified ERA file over the area specified by the dem rsc file. ''' def __init__(self,gribfile,demfile,grib='ECMWF',humidity='Q',demtype=np.int16,demfmt='HGT',verb=False, Del='comb'): ''' Initiates the data structure for atmos corrections in geocoded domain. Args: * gribfile (str) : Path to the weather model file. * demfile (str) : Path to the SIM_nRLKS.hgt file. Kwargs: * grib (str) : Can be ECMWF, ERA, NARR or MERRA depending on the model. * humidity (str) : Can be 'Q' or 'R' depdending on specific or relative humidity. * demtype (dtype) : Data type of the DEM File. * demfmt (str) : Can be 'HGT' or 'RMG' depending on type of file. If RMG, altitude is assumed to be second channel. * Del (str) : Can be 'comb', 'dry' or 'wet'. Specifies which delay you want (both, dry or wet) .. note :: The DEMFile is only used to set up the geographic boundaries of the problem and the size of the output matrix during init. The GRIBfile is read and the 3D interpolation function is prepared. ''' sanity_check(grib) grib = grib.upper() humidity = humidity.upper() demfmt = demfmt.upper() assert grib in ('ERA','NARR','ECMWF','MERRA'), 'PyAPS: Undefined grib file source.' self.grib = grib '''String with the name of the weather model.''' assert humidity in ('Q','R'), 'PyAPS: Undefined humidity.' self.humidity = humidity if self.grib in ('NARR','MERRA'): assert self.humidity in ('Q'), 'PyAPS: Relative humidity not provided by NARR/MERRA.' assert os.path.isfile(gribfile), 'PyAPS: GRIB File does not exist.' self.gfile = gribfile assert os.path.isfile(demfile), 'PyAPS: DEM file does not exist.' self.hfile = demfile assert demfmt in ('RMG','HGT','VAR'), 'PyAPS: DEM Format can be RMG or HGT, VAR' self.fmt = demfmt self.demtype = demtype self.dict = processor.initconst() self.bufspc = 1.2 ####Hardcoded. if grib in ('ERA','ECMWF'): self.hgtscale = ((self.dict['maxAlt']-self.dict['minAlt'])/self.dict['nhgt'])/0.703 elif grib in ('NARR'): self.hgtscale = ((self.dict['maxAlt']-self.dict['minAlt'])/self.dict['nhgt'])/0.3 elif grib in ('MERRA'): self.hgtscale = ((self.dict['maxAlt']-self.dict['minAlt'])/self.dict['nhgt'])/0.5 [lon,lat,nx,ny] = utils.geo_rsc(self.hfile,verbose=verb) ######Problems near the international dateline if grib in ('MERRA'): #MERRA's longtidue range -180~180. #This is to correct the value from utils.geo_rsc() function. Cunren, OCT-2015 lon[lon > 180] -= 360.0 #if(lon[0] > 180): # lon[0] = lon[0] - 360.0 #if(lon[1] > 180): # lon[1] = lon[1] - 360.0 self.minlon = lon.min()-self.bufspc self.maxlon = lon.max()+self.bufspc self.minlat = lat.min()-self.bufspc self.maxlat = lat.max()+self.bufspc self.nx = nx self.ny = ny else: lon[lon < 0] += 360.0 self.minlon = lon.min()-self.bufspc self.maxlon = lon.max()+self.bufspc self.minlat = lat.min()-self.bufspc self.maxlat = lat.max()+self.bufspc self.nx = nx self.ny = ny if self.grib in ('ERA'): [lvls,latlist,lonlist,gph,tmp,vpr] = era.get_era(self.gfile, self.minlat, self.maxlat, self.minlon, self.maxlon, self.dict, humidity=self.humidity, verbose=verb) elif self.grib in ('ECMWF'): [lvls,latlist,lonlist,gph,tmp,vpr] = era.get_ecmwf(self.gfile, self.minlat, self.maxlat, self.minlon, self.maxlon, self.dict, humidity=self.humidity, verbose=verb) elif self.grib in ('NARR'): [lvls,latlist,lonlist,gph,tmp,vpr] = narr.get_narr(self.gfile, self.minlat, self.maxlat, self.minlon, self.maxlon, self.dict, verbose=verb) elif self.grib in ('MERRA'): [lvls,latlist,lonlist,gph,tmp,vpr] = merra.get_merra(self.gfile, self.minlat, self.maxlat, self.minlon, self.maxlon, self.dict, verbose=verb) lonlist[lonlist < 0.] += 360.0 hgt = np.linspace(self.dict['minAlt'],self.dict['maxAlt'],self.dict['nhgt']) [Pi,Ti,Vi] = processor.intP2H(lvls,hgt,gph,tmp,vpr,self.dict,verbose=verb) [DDry,DWet] = processor.PTV2del(Pi,Ti,Vi,hgt,self.dict,verbose=verb) if Del in ('comb','Comb'): Delfn=DDry+DWet elif Del in ('dry','Dry'): Delfn=DDry elif Del in ('wet','Wet'): Delfn=DWet else: print('Unrecognized delay type') sys.exit(1) if self.grib in ('NARR'): [Delfn,latlist,lonlist] = narr.intdel(hgt,latlist,lonlist,Delfn) self.Delfn = Delfn self.lonlist = lonlist self.latlist = latlist self.Pi = Pi self.Vi = Vi self.Ti = Ti self.hgt = hgt self.verb = verb def merisfactor(self,dataobj,inc=0.0,wvl=4*np.pi): ''' Write pi-factor from Li et al 2012 to a matrix / HDF5 object or a file directly. Args: * dataobj (str or HDF5 or np.array): Final output. If str, output is written to file. Kwargs: * inc (np.float): Incidence angle in degrees. Default is vertical. * wvl (np.float): Wavelength in meters. Default output results in delay in meters. .. note:: If dataobj is string, output is written to the file. If np.array or HDF5 object, it should be of size (ny,nx). ''' minAltp = self.dict['minAltP'] # Incidence cinc = np.cos(inc*np.pi/180.0) # Compute the two integrals WonT = self.Vi/self.Ti WonT2 = WonT/self.Ti S1 = intg.cumtrapz(WonT,x=self.hgt,axis=-1) val = 2*S1[:,-1]-S1[:,-2] val = val[:,None] S1 = np.concatenate((S1,val),axis=-1) del WonT S2 = intg.cumtrapz(WonT2,x=self.hgt,axis=-1) val = 2*S2[:,-1]-S2[:,-2] val = val[:,None] S2 = np.concatenate((S2,val),axis=-1) del WonT2 Tm = S1/S2 self.Tm = Tm # Reading in the DEM if self.verb: print('PROGRESS: READING DEM') fin = open(self.hfile,'rb'); if self.fmt in ('HGT'): dem = np.fromfile(file=fin,dtype=self.demtype,count=self.nx*self.ny).reshape(self.ny,self.nx) elif self.fmt in ('RMG'): dem = np.fromfile(file=fin,dtype=self.demtype,count=2*self.nx*self.ny).reshape(self.ny,2*self.nx) dem = dem[:,self.nx:] dem = np.round(dem).astype(np.int) fin.close() # check output, and open file if necessary outFile = isinstance(dataobj,str) if outFile: fout = open(dataobj,'wb') dout = np.zeros((self.ny,self.nx)) else: assert ((dataobj.shape[0]==self.ny) & (dataobj.shape[1]==self.nx)), 'PyAPS: Not a valid data object.' dout = dataobj # Create the lon/lat arrays laty = np.linspace(self.maxlat-self.bufspc,self.minlat+self.bufspc,self.ny) lonx = np.linspace(self.minlon+self.bufspc,self.maxlon-self.bufspc,self.nx) # Create the 1d interpolator if self.verb: print('PROGRESS: FINE INTERPOLATION OF HEIGHT LEVELS') intp_1d = si.interp1d(self.hgt,Tm,kind='cubic',axis=1) # Interpolate the Tm variable every meter dem[dem < minAltp] = minAltp minH = dem.min() maxH = dem.max()+1 kh = np.arange(minH,maxH) Tm_1m = intp_1d(kh) self.alti = kh # Reshape Tm Lonu = np.unique(self.lonlist) Latu = np.unique(self.latlist) nLon = len(Lonu) nLat = len(Latu) Tm_1m = np.reshape(Tm_1m.T,(len(kh),nLat,nLon)) self.Tm_1m = Tm_1m # Create the cube interpolator for the bilinear method if self.verb: print('PROGRESS: CREATE THE BILINEAR INTERPOLATION FUNCTION') bilicube = processor.Bilinear2DInterpolator(Lonu,Latu,Tm_1m,cube=True) # Get the values from the dictionnary k1 = self.dict['k1'] k2 = self.dict['k2'] k3 = self.dict['k3'] mmO = self.dict['mmO'] mmH = self.dict['mmH'] mma = self.dict['mma'] w = (2*mmH + mmO)/mma Rv = self.dict['Rv'] Rho = self.dict['Rho'] # Loop on the lines if self.verb: toto = utils.ProgressBar(maxValue=self.ny) print('PROGRESS: MAPPING THE DELAY') for m in range(self.ny): if self.verb: toto.update(m,every=5) # Get Lon/Lat loni = lonx lati = laty[m]*np.ones((loni.shape)) # Make the bilinear interpolation D = dem[m,:] - minH val = bilicube(loni,lati,D) val = 0.000001 * Rho * Rv * ( k3/val + k2 - w*k1) * np.pi*4.0/(cinc*wvl) if outFile: resy = val.astype(np.float32) resy.tofile(fout) else: dataobj[m,:] = val if self.verb: toto.close() # Close if outfile if outFile: fout.close() ###########End of PyAPS_geo class###############################