From b82574216b8a0bb934eabd81629e11910db80d2a Mon Sep 17 00:00:00 2001 From: Ryan Torn Date: Fri, 29 Jul 2022 16:54:57 +0000 Subject: [PATCH] removing files in another repo --- ecmwf_teton.py | 381 ---------------------------------------------- grid_calc.f90 | 145 ------------------ tigge_download.py | 269 -------------------------------- 3 files changed, 795 deletions(-) delete mode 100644 ecmwf_teton.py delete mode 100644 grid_calc.f90 delete mode 100644 tigge_download.py diff --git a/ecmwf_teton.py b/ecmwf_teton.py deleted file mode 100644 index 6dd8f3a..0000000 --- a/ecmwf_teton.py +++ /dev/null @@ -1,381 +0,0 @@ -import os -import time -import shutil -import sys -import cfgrib -import datetime as dt -import glob -import gzip -import urllib -import pandas as pd -import numpy as np -import xarray as xr - -def stage_grib_files(datea, config): - ''' - This is a generic class for copying or linking grib file data from a specific location - to a directory where calculations are performed. No matter where these calculations are - carried out, this routine must exist. - - This particular instance is for ECMWF data on the machine teton at UAlbany. In this - case, the grib files are linked to files that are located in another directory on that - machine. - - Attributes: - datea (string): The initialization time of the forecast (yyyymmddhh) - config (dict): The dictionary with configuration information - ''' - - freq = config.get('fcst_hour_int', 12) - fmax = config.get('fcst_hour_max', 120) - - # Make the work directory if it does not exist - if not os.path.isdir(config['work_dir']): - try: - os.makedirs(config['work_dir']) - except OSError as e: - raise e - - init = dt.datetime.strptime(datea, '%Y%m%d%H') - init_s = init.strftime("%m%d%H%M") - - # Loop over all forecast times, link to the source file - for fhr in range(0, int(fmax)+int(freq), int(freq)): - - datef = init + dt.timedelta(hours=fhr) - datef_s = datef.strftime("%m%d%H%M") - - grib_file = 'E1E{0}{1}1'.format(str(init_s), str(datef_s)) - infile = '{0}/{1}'.format(config['model_dir'],grib_file) - outfile = '{0}/{1}'.format(config['work_dir'],grib_file) - - # Only try to copy if the file is not there - if ( not os.path.isfile(outfile) ): - - # Wait for the source file to be present - while not os.path.exists(infile): - time.sleep(20.1) - - # Wait for the file to be finished being copied - while ( (time.time() - os.path.getmtime(infile)) < 60 ): - time.sleep(10) - - try: # Try to link from the source to the work directory - os.symlink(infile, outfile) - except Exception as err: - raise err - - -def stage_atcf_files(datea, bbnnyyyy, config): - ''' - This is a generic routine for copying or linking ATCF file data from a specific location - to a directory where calculations are performed. No matter where these calculations are - carried out, this routine must exist. - - The result is a set of ATCF files in the work directory of the format atcf_NN.dat, - where NN is the ensemble member number. - - This particular instance is for ECMWF data on the machine teton at UAlbany. In this - case, all ATCF data for a particular storm is in one file, so the code waits for this - initialization time to exist, then uses sed to get the lines attributed to each - ensemble member and places that data in a seperate file. - - Attributes: - datea (string): The initialization time of the forecast (yyyymmddhh) - bbnnyyyy (string): TC Identification, where bb is the basin, nn is the number, yyyy is year - config (dict): The dictionary with configuration information - ''' - - src = '{0}/a{1}.dat'.format(config['atcf_dir'],bbnnyyyy) - nens = int(config['num_ens']) - - # Wait for the source file to be present - while not os.path.exists(src): - time.sleep(20.5) - - # Wait for the ensemble ATCF information to be placed in the file - while ( len(os.popen('sed -ne /{0}/p {1} | sed -ne /EE/p'.format(datea,src)).read()) == 0 ): - time.sleep(20.7) - - # Wait for the file to be finished being copied - while ( (time.time() - os.path.getmtime(src)) < 60 ): - time.sleep(10) - - for n in range(nens + 1): - - nn = '%0.2i' % n - file_name = '{0}/atcf_{1}.dat'.format(config['work_dir'],nn) - - # If the specific member's ATCF file does not exist, copy from the source file with sed. - if not os.path.isfile(file_name): - - fo = open(file_name,"w") - fo.write(os.popen('sed -ne /{0}/p {1} | sed -ne /EE{2}/p'.format(datea,src,nn)).read()) - fo.close() - - -def stage_best_file(bbnnyyyy, config): - ''' - This is a generic routine for copying or linking ATCF best track file data from a specific location - to a directory where calculations are performed. No matter where these calculations are - carried out, this routine must exist. - - The result is a single best track file in the work directory of the format bbnnyyyy.dat, - where bb is the basin, nn is the TC number. - - This particular instance copies data from the NHC data server. - - Attributes: - bbnnyyyy (string): TC Identification, where bb is the basin, nn is the number, yyyy is year - config (dict): The dictionary with configuration information - ''' - - try: # Look for the data in the real-time directory - - filei = urllib.request.urlopen('{0}/b{1}.dat'.format(config.get('best_dir','https://ftp.nhc.noaa.gov/atcf/btk'),bbnnyyyy)) - fileo = open('{0}/b{1}.dat'.format(config['work_dir'],bbnnyyyy), 'wb') - fileo.write(filei.read()) - filei.close() - fileo.close() - - except: # If the file is not present in the real-time directory, look in the archive. - - src = '{0}/{1}/b{2}.dat.gz'.format(config.get('best_dir_alt','https://ftp.nhc.noaa.gov/atcf/archive'),bbnnyyyy[4:8],bbnnyyyy) - - # Unzip the file from the NHC server, write the file to the work directory - gzfile = gzip.GzipFile(fileobj=urllib.request.urlopen(src)) - uzfile = open('{0}/b{1}.dat'.format(config['work_dir'],bbnnyyyy), 'wb') - uzfile.write(gzfile.read()) - gzfile.close() - uzfile.close() - - -class ReadGribFiles: - ''' - This is a generic class that is used to read specific forecast fields from a grib file at a specific - forecast hour. This class is a generic class, but differs depending on the grib file format and fields. - The init routine creates the class and constructs a field dictionary. - - This particular instance is for ECMWF data on the UAlbany teton machine. Each ECMWF file contains all - members for a specific forecast time. - - Attributes: - datea (string): The initialization time of the forecast (yyyymmddhh) - fhr (int): Forecast hour - config (dict): The dictionary with configuration information - ''' - - def __init__(self, datea, fhr, config): - - init = dt.datetime.strptime(datea, '%Y%m%d%H') - init_s = init.strftime("%m%d%H%M") - datef = init + dt.timedelta(hours=fhr) - datef_s = datef.strftime("%m%d%H%M") - - # Construct the grib file dictionary for a particular forecast hour - file_name = os.path.join(config['work_dir'], "E1E{0}{1}1".format(str(init_s), str(datef_s))) - try: - ds = cfgrib.open_datasets(file_name) - self.grib_dict = {} - for d in ds: - for tt in d: - if 'number' in d[tt].dims: - self.grib_dict.update({'{0}_pf'.format(tt): d[tt]}) - else: - self.grib_dict.update({'{0}_cf'.format(tt): d[tt]}) - - except IOError as exc: - raise RuntimeError('Failed to open {0}'.format(file_name)) from exc - - # This is a dictionary that maps from generic variable names to the name of variable in file - self.var_dict = {'zonal_wind': 'u', \ - 'meridional_wind': 'v', \ - 'zonal_wind_10m': 'u10', \ - 'meridional_wind_10m': 'v10', \ - 'geopotential_height': 'gh', \ - 'temperature': 't', \ - 'relative_humidity': 'r', \ - 'specific_humidity': 'q', \ - 'sea_level_pressure': 'msl', \ - 'precipitation': 'tp'} - - for key in self.grib_dict: - if np.max(self.grib_dict[key].coords['longitude']) > 180: - self.grib_dict[key].coords['longitude'] = (self.grib_dict[key].coords['longitude'] + 180) % 360 - 180 - - if config.get('flip_lon','False') == 'True': - for key in self.grib_dict: - self.grib_dict[key].coords['longitude'] = (self.grib_dict[key].coords['longitude'] + 360.) % 360. - self.grib_dict[key] = self.grib_dict[key].sortby('longitude') - - if '{0}_cf'.format(self.var_dict['specific_humidity']) in self.grib_dict: - self.has_specific_humidity = True - else: - self.has_specific_humidity = False - - self.has_total_precip = True - - self.nens = int(self.grib_dict['gh_pf'].attrs['GRIB_totalNumber']) - - - def set_var_bounds(self, varname, vdict): - ''' - This is a generic routine that is used to determine the appropriate starting and ending latitude, - longitude, and if appropriate pressure levels to extract from the grib files. This routine takes into - account the order of the coordinate variables. If the bounds are not specified, the routine uses all - coordinate indicies. The resulting bounds are added back into the dictionary object and returned for - use in other routines within the class. - - Attributes: - varname (string): The name of the variable that will be extracted from file - vdict (dict): The dictionary object with variable information - ''' - - vname = '{0}_cf'.format(self.var_dict[varname]) - - if 'latitude' in vdict: - - # See if the latitude values are ordered reversed - if float(self.grib_dict[vname].attrs['GRIB_latitudeOfFirstGridPointInDegrees']) > \ - float(self.grib_dict[vname].attrs['GRIB_latitudeOfLastGridPointInDegrees']): - vdict['lat_start'] = int(vdict['latitude'][1]) - vdict['lat_end'] = int(vdict['latitude'][0]) - else: - vdict['lat_start'] = int(vdict['latitude'][0]) - vdict['lat_end'] = int(vdict['latitude'][1]) - - else: - - # Get all latitude values - latvec = list(self.grib_dict[vname].latitude.data) - vdict['lat_start'] = latvec[0] - vdict['lat_end'] = latvec[-1] - - if 'longitude' in vdict: - - # Take longitude from input values - vdict['lon_start'] = int(vdict['longitude'][0]) - vdict['lon_end'] = int(vdict['longitude'][1]) - - else: - - # Use all longitude values - lonvec = list(self.grib_dict[vname].longitude.data) - vdict['lon_start'] = lonvec[0] - vdict['lon_end'] = lonvec[-1] - - if 'isobaricInhPa' in vdict: - - # See of pressure level values are reversed - if self.grib_dict[vname].isobaricInhPa[0] > self.grib_dict[vname].isobaricInhPa[1]: - vdict['pres_start'] = int(vdict['isobaricInhPa'][1]) - vdict['pres_end'] = int(vdict['isobaricInhPa'][0]) - else: - vdict['pres_start'] = int(vdict['isobaricInhPa'][0]) - vdict['pres_end'] = int(vdict['isobaricInhPa'][1]) - - return vdict - - - def create_ens_array(self, varname, nens, vdict): - ''' - This is a generic routine that is used to create an xarray object that contains all ensemble - members for a particular field, with the proper units and coordinate variables. The resulting - array has dimensions (ensemble, latitude, longitude). The routine - takes into account the order of the lat/lon arrays - - Attributes: - varname (string): The name of the variable that will be extracted from file - nens (int): number of ensemble members - vdict (dict): The dictionary object with variable information - ''' - - vname = '{0}_cf'.format(self.var_dict[varname]) - - # Create attributes based on what is in the file - attrlist = {} - if 'description' in vdict: - attrlist['description'] = vdict['description'] - if 'units' in vdict: - attrlist['units'] = vdict['units'] - if '_FillValue' in vdict: - attrlist['_FillValue'] = vdict['_FillValue'] - - # Create an empty xarray that can be used to copy data into - lonvec = list(self.grib_dict[vname].sel(longitude=slice(vdict['lon_start'], vdict['lon_end'])).longitude.data) - latvec = list(self.grib_dict[vname].sel(latitude=slice(vdict['lat_start'], vdict['lat_end'])).latitude.data) - ensarr = xr.DataArray(name='ensemble_data', data=np.zeros([nens, len(latvec), len(lonvec)]), - dims=['ensemble', 'latitude', 'longitude'], attrs=attrlist, - coords={'ensemble': [i for i in range(nens)], 'latitude': latvec, 'longitude': lonvec}) - - return(ensarr) - - - def read_pressure_levels(self, varname): - - vname = '{0}_cf'.format(self.var_dict[varname]) - return self.grib_dict[vname].isobaricInhPa[:] - - - # Function to read a single ensemble member's forecast field - def read_grib_field(self, varname, member, vdict): - ''' - This is a generic routine that is used to read a forecast field from a single ensemble member - based on the information contained within the dictionary vdict. - - Attributes: - varname (string): Name of the variable that will be extracted from file - member (int): ensemble member to read from file - vdict (dict): The dictionary object with variable information - ''' - - # Read a single pressure level of data, if this is a variable that has pressure levels - if 'isobaricInhPa' in vdict: - - if member == 0: # Control member - vname = '{0}_cf'.format(self.var_dict[varname]) - vout = self.grib_dict[vname].sel(latitude=slice(vdict['lat_start'], vdict['lat_end']), \ - longitude=slice(vdict['lon_start'], vdict['lon_end']), \ - isobaricInhPa=slice(vdict['pres_start'], vdict['pres_end'])) - - else: - vname = '{0}_pf'.format(self.var_dict[varname]) - vout = self.grib_dict[vname].sel(number=member, \ - latitude=slice(vdict['lat_start'], vdict['lat_end']), \ - longitude=slice(vdict['lon_start'], vdict['lon_end']), \ - isobaricInhPa=slice(vdict['pres_start'], vdict['pres_end'])) - - # Read the only level if it is a single level variable - else: - - if member == 0: # Control member - vname = '{0}_cf'.format(self.var_dict[varname]) - vout = self.grib_dict[vname].sel(latitude=slice(vdict['lat_start'], vdict['lat_end']), \ - longitude=slice(vdict['lon_start'], vdict['lon_end'])) - - else: - - vname = '{0}_pf'.format(self.var_dict[varname]) - vout = self.grib_dict[vname].sel(number=member, \ - latitude=slice(vdict['lat_start'], vdict['lat_end']), \ - longitude=slice(vdict['lon_start'], vdict['lon_end'])) - - return(vout) - - - def close_files(self): - - del self.grib_dict - - -if __name__ == '__main__': - - src1 = "/Users/parthpatwari/RA_Atmospheric_Science/Old_Code/atcf_data" - grib_src = "/Users/parthpatwari/RA_Atmospheric_Science/GRIB_files" - dest1 = "/Users/parthpatwari/RA_Atmospheric_Science/New_Code/atcf_data" - atcf_src = "/Users/parthpatwari/RA_Atmospheric_Science/Old_Code/atcf_data" - # c1 = CopyFiles(src, dest) - # if c1.checkandcreatedir(): - # c1.copy_filestowork() - g1 = ReadGribFiles(grib_src, '2019082900', 180) - a1 = Readatcfdata(atcf_src) diff --git a/grid_calc.f90 b/grid_calc.f90 deleted file mode 100644 index 1586e1d..0000000 --- a/grid_calc.f90 +++ /dev/null @@ -1,145 +0,0 @@ -! python -m numpy.f2py -c grid_calc.f90 -m grid_calc - -subroutine calc_circ(vort, circrad, dx, nx, ny, circ) - - implicit none - - integer, intent(in) :: nx, ny - real, intent(in) :: vort(ny,nx), circrad, dx - real, intent(out) :: circ(ny,nx) - - integer :: cint, i, ii, j, jj - real :: csum, asum - - cint = nint(circrad / dx) - - do ii = 1, nx - do jj = 1, ny - - csum = 0.0 - asum = 0.0 - do i = max(ii-cint, 1), min(ii+cint, nx) - do j = max(jj-cint, 1), min(jj+cint, ny) - if ( sqrt(real(ii-i)**2 + real(jj-j)**2)*dx <= circrad ) then - csum = csum + vort(j,i)*dx*dx - asum = asum + dx*dx - end if - end do - end do - circ(jj,ii) = csum / asum - - end do - end do - - return -end subroutine calc_circ - - -subroutine calc_circ_llgrid(vort, circrad, lat, lon, global, nx, ny, circ) - - implicit none - - real, parameter :: pid = atan(1.0) / 45.0 - real, parameter :: Rearth = 6378.388 - - integer, intent(in) :: nx, ny - logical, intent(in) :: global - real, intent(in) :: vort(ny,nx), lat(ny), lon(nx), circrad - real, intent(out) :: circ(ny,nx) - real, allocatable :: dgrid(:,:) - - integer :: xint, yint, i, ii, iii, j, jj, jjj, k, i1, i2, j1, j2, ngy - real :: csum, asum, abox(ny,nx), dlon, dlat, gc_dist - - dlon = abs(lon(2)-lon(1)) - dlat = abs(lat(2)-lat(1)) - yint = ceiling(circrad / (Rearth*pid*dlat)) - - ! Compute the area of a latitutde/longitude box - do i = 1, nx - do j = 1, ny - abox(j,i) = abs(dlon*pid*cos(pid*lat(j))*dlat*pid) - end do - end do - - do jj = 1, ny - - xint = min(ceiling(abs(circrad / (dlon*pid*Rearth*cos(pid*lat(jj))))),nx/2) - - j1 = max(jj-yint,1) - j2 = min(jj+yint,ny) - ngy = j2-j1+1 - - ! Figure out all of the points that are within radius circrad - allocate(dgrid(ngy,2*xint+1)) - do i = 1, 2*xint+1 - do j = 1, ngy - if ( gc_dist(lat(j+j1-1),lon(i),lat(jj),lon(xint+1)) <= circrad ) then - dgrid(j,i) = 1.0 - else - dgrid(j,i) = 0.0 - end if - end do - end do - - do ii = 1, nx - - csum = 0.0 - asum = 0.0 - - if ( global ) then - i1 = ii - xint - i2 = ii + xint - else - i1 = max(ii-xint, 1) - i2 = min(ii+xint, nx) - end if - - ! Loop over all points near a specific center point - do i = 1, i2-i1+1 - - iii = i + i1 - 1 - if ( iii < 1 .and. global ) iii = iii + nx - if ( iii > nx .and. global ) iii = iii - nx - - ! Compute the area-average quantity and total area - do j = 1, ngy - jjj = j + j1 - 1 - csum = csum + vort(jjj,iii)*abox(jjj,iii)*dgrid(j,i) - asum = asum + abox(jjj,iii)*dgrid(j,i) - end do - - end do - - ! Divide by area to get area-average vorticity - if ( asum > 0.0 ) then - circ(jj,ii) = csum / asum - end if - - end do - - deallocate(dgrid) - - end do - - return -end subroutine calc_circ_llgrid - - -function gc_dist(lat1,lon1,lat2,lon2) - -implicit none - -real, parameter :: pid = atan(1.0) / 45.0 -real, parameter :: Rearth = 6378.388 - -real, intent(in) :: lat1, lon1, lat2, lon2 - -real :: gc_dist - -gc_dist = Rearth * acos( min(sin(lat1*pid) * sin(lat2*pid) + & - cos(lat1*pid) * cos(lat2*pid) * & - cos((lon2-lon1)*pid), 1.0) ) - -return -end function gc_dist diff --git a/tigge_download.py b/tigge_download.py deleted file mode 100644 index c09156c..0000000 --- a/tigge_download.py +++ /dev/null @@ -1,269 +0,0 @@ -from ecmwfapi import ECMWFDataServer -import sys, os -import argparse -import datetime as dt -import configparser -import numpy as np - -''' -Program that retrieves forecast fields from a single initialization time from the -TIGGE database. - -From command line: - -python run_AR_sens.py -init yyyymmddhh --param paramfile - -where: - - -init is the initialization date in yyyymmddhh format - -storm is TC name, where XXXXXXXNNB, where XXXXXXXX is the name, NN is the number, B is the basin (optional) - -param is the parameter file path (optional, otherwise goes to default values in example.parm) -''' - -# Read the initialization time and storm from the command line -exp_parser = argparse.ArgumentParser() -exp_parser.add_argument('--init', action='store', type=str, required=True) -exp_parser.add_argument('--storm', action='store', type=str) -exp_parser.add_argument('--param', action='store', type=str) - -args = exp_parser.parse_args() - -yyyymmddhh = args.init - -# Read the parameter file -if args.param: - paramfile = args.param -else: - paramfile = 'example.parm' - -confin = configparser.ConfigParser() -confin.read(paramfile) - -config = {} -config.update(confin['model']) -config.update(confin['locations']) - -# Modify work and output directory for specific case/time -if args.storm: - storm = args.storm - config['work_dir'] = '{0}/{1}.{2}'.format(config['work_dir'],yyyymmddhh,storm) -else: - config['work_dir'] = '{0}/{1}'.format(config['work_dir'],yyyymmddhh) - -# Create appropriate directories -if not os.path.isdir(config['work_dir']): - try: - os.makedirs(config['work_dir']) - except OSError as e: - raise e - -# Set the TIGGE read command values -if not config['tigge_forecast_time']: - raise "tigge_forecast_time is missing" - -daystr = yyyymmddhh[0:4] + '-' + yyyymmddhh[4:6] + '-' + yyyymmddhh[6:8] -hhstr = yyyymmddhh[8:10] + ":00:00" - -if config['model_src'] == 'ECMWF': - model_str = "ecmf" - model_num = "1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20/21/22/23/24/25/26/27/28/29/30/31/32/33/34/35/36/37/38/39/40/41/42/43/44/45/46/47/48/49/50" -elif config['model_src'] == 'GEFS': - model_str = "kwbc" - model_num = "1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20/21/22/23/24/25/26/27/28/29/30" -# model_num = "1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20" - -atmres = config.get('tigge_forecast_grid_space','1.0/1.0') -sfcres = "0.25/0.25" - -os.chdir(config['work_dir']) - -server = ECMWFDataServer() - -# Retrieve pressure level data -server.retrieve({ - 'origin' : model_str, - 'levelist' : "200/250/300/500/700/850/925/1000", - 'levtype' : "pl", - 'expver' : "prod", - 'parameter' : "130/131/132/133/156", - 'number' : model_num, - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : atmres, - 'time' : hhstr, - 'type' : "pf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_pl_pf.grib" -}) - -server.retrieve({ - 'origin' : model_str, - 'levelist' : "200/250/300/500/700/850/925/1000", - 'levtype' : "pl", - 'expver' : "prod", - 'parameter' : "130/131/132/133/156", - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : atmres, - 'time' : hhstr, - 'type' : "cf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_pl_cf.grib" -}) - -''' -# Retrieve PV level data -server.retrieve({ - 'origin' : model_str, - 'levtype' : "pv", - 'levelist' : "2", - 'expver' : "prod", - 'parameter' : "3/131/132", - 'number' : model_num, - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : atmres, - 'time' : hhstr, - 'type' : "pf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_pv_pf.grib" -}) - -server.retrieve({ - 'origin' : model_str, - 'levtype' : "pv", - 'levelist' : "2", - 'expver' : "prod", - 'parameter' : "3/131/132", - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : atmres, - 'time' : hhstr, - 'type' : "cf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_pv_cf.grib" -}) - -''' - -# Retrieve surface data -server.retrieve({ - 'origin' : model_str, - 'levtype' : "sfc", - 'expver' : "prod", - 'parameter' : "59/136/151/167/168", - 'number' : model_num, - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : atmres, - 'time' : hhstr, - 'type' : "pf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_sfc_pf.grib" -}) - -server.retrieve({ - 'origin' : model_str, - 'levtype' : "sfc", - 'expver' : "prod", - 'parameter' : "59/136/151/167/168", - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : atmres, - 'time' : hhstr, - 'type' : "cf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_sfc_cf.grib" -}) - -# Retrieve higher-resolution surface fields -server.retrieve({ - 'origin' : model_str, - 'levtype' : "sfc", - 'expver' : "prod", - 'parameter' : "165/166/228", - 'number' : model_num, - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : sfcres, - 'time' : hhstr, - 'type' : "pf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_hrsfc_pf.grib" -}) - -server.retrieve({ - 'origin' : model_str, - 'levtype' : "sfc", - 'expver' : "prod", - 'parameter' : "165/166/228", - 'dataset' : "tigge", - 'step' : config['tigge_forecast_time'], - 'grid' : sfcres, - 'time' : hhstr, - 'type' : "cf", - 'date' : daystr, - 'class' : "ti", - 'target' : "tigge_output_hrsfc_cf.grib" -}) - -ftype = ("pl", "sfc", "hrsfc") - -for i in range(len(ftype)): - os.system("wgrib2 -s tigge_output_{0}_pf.grib >& grib_{0}_pf.out".format(ftype[i])) - os.system("wgrib2 -s tigge_output_{0}_cf.grib >& grib_{0}_cf.out".format(ftype[i])) - -hlist = config['tigge_forecast_time'].split("/") - -init = dt.datetime.strptime(yyyymmddhh, '%Y%m%d%H') -init_s = init.strftime("%m%d%H%M") - -# Loop over all times, create one file per time -for t in range(len(hlist)): - - datef = init + dt.timedelta(hours=int(hlist[t])) - datef_s = datef.strftime("%m%d%H%M") - - if int(hlist[t]) == 0: - timestr = ":anl:" - else: - timestr = ":{0} hour fcst:".format(hlist[t]) - hhh = '%0.3i' % int(hlist[t]) - - if np.remainder(int(hlist[t]),24) == 0: - pcpstr = ":0-{0} day acc fcst:".format(round(float(hlist[t]) / 24.)) - else: - pcpstr = ":0-{0} hour acc fcst:".format(hlist[t]) - - if os.path.isfile("f{0}_fields.grb".format(hhh)): - os.remove("f{0}_fields.grb".format(hhh)) - - print(hlist[t],hhh,timestr,pcpstr) - - gribout = 'E1E{0}{1}1'.format(str(init_s), str(datef_s)) - - for i in range(len(ftype)): - os.system("cat grib_{0}_cf.out | grep \"{1}\" | wgrib2 -fix_ncep -i -append tigge_output_{0}_cf.grib \ - -grib {2} >& /dev/null".format(ftype[i],timestr,gribout)) - os.system("cat grib_{0}_pf.out | grep \"{1}\" | wgrib2 -fix_ncep -i -append tigge_output_{0}_pf.grib \ - -grib {2} >& /dev/null".format(ftype[i],timestr,gribout)) - - if int(hlist[t]) > 0: - os.system("cat grib_hrsfc_cf.out | grep \"{0}\" | wgrib2 -fix_ncep -i -append tigge_output_hrsfc_cf.grib \ - -grib {1} >& /dev/null".format(pcpstr,gribout)) - os.system("cat grib_hrsfc_pf.out | grep \"{0}\" | wgrib2 -fix_ncep -i -append tigge_output_hrsfc_pf.grib \ - -grib {1} >& /dev/null".format(pcpstr,gribout)) - -for i in range(len(ftype)): - os.remove("tigge_output_{0}_pf.grib".format(ftype[i])) - os.remove("tigge_output_{0}_cf.grib".format(ftype[i])) - os.remove("grib_{0}_pf.out".format(ftype[i])) - os.remove("grib_{0}_cf.out".format(ftype[i])) -