Skip to content
Snippets Groups Projects
Select Git revision
  • 68fc2efec0023076d7626b256bbc64cc8d3627a3
  • consistent_config default protected
2 results

generate_nature.py

Blame
  • create_obsseq.py 11.67 KiB
    """The functions in here create obs_seq.in files.
    These are the template files defining obs locations and metadata
    according to which observations are generated and subsequently assimilated.
    """
    import os, sys, warnings
    import numpy as np
    import datetime as dt
    from pysolar.solar import get_altitude, get_azimuth
    
    from config.cfg import exp
    from config.cluster import cluster
    from dartwrf.obs import calculate_obs_locations as col
    
    def obskind_read():
        """Read dictionary of observation types + ID numbers ("kind") 
        from DART f90 script
        """
    
        definitionfile = cluster.scriptsdir+'/../config/obs_kind_mod.f90'
        with open(definitionfile, 'r') as f:
            kind_def_f = f.readlines()
    
        obskind_nrs = {}
        for i, line in enumerate(kind_def_f):
            if 'Integer definitions for DART OBS TYPES' in line:
                # data starts below this line
                i_start = i
                break
        for line in kind_def_f[i_start+1:]:
            if 'MAX_DEFINED_TYPES_OF_OBS' in line:
                # end of data
                break
            if '::' in line:
                # a line looks like this
                # integer, parameter, public ::       MSG_4_SEVIRI_TB =   261
                data = line.split('::')[-1].split('=')
                kind_str = data[0].strip()
                kind_nr = int(data[1].strip())
                obskind_nrs[kind_str] = kind_nr
        return obskind_nrs
    
    
    #####################
    # Global variables
    obs_kind_nrs = obskind_read()  # DART internal indices
    
    # position on earth for RTTOV ray geometry
    lat0 = 45.
    lon0 = 0.
    sat_az = "180.0"
    sat_zen = "45.0"
    
    
    def degr_to_rad(degr):
        """Convert to DART convention (radians)
        2*pi = 360 degrees
        
        Args:
            degr (float) : degrees east of Greenwich
            
        Returns 
            float
        """
        if degr < 0:
            degr += 360
        return degr/360*2*np.pi
    
    
    def round_to_day(dtobj):
        """Overwrite hours, minutes, seconds to 0
        Args:
            dtobj (dt.datetime)
        """
        return dtobj.replace(second=0, minute=0, hour=0)
    
    
    def add_timezone_UTC(t):
        return dt.datetime(t.year, t.month, t.day, t.hour, t.minute, tzinfo=dt.timezone.utc)
    
    
    def get_dart_date(time_dt):
        """Convert datetime.datetime into DART time format
        
        Assumes that input is UTC!
        
        Returns
            str, str
        """
        days_since_010120 = 145731
        ref_day = dt.datetime(2000, 1, 1, tzinfo=dt.timezone.utc)
        dart_date_day = str((time_dt - ref_day).days + days_since_010120)
        secs_thatday = str(int((time_dt - round_to_day(time_dt)).total_seconds()))
        return dart_date_day, secs_thatday
    
    
    def write_tuple_to_pickle(fpath_out, tuple):
        import pickle
        os.makedirs(os.path.dirname(fpath_out), exist_ok=True)
        with open(fpath_out, 'wb') as f:
            pickle.dump(tuple, f)
        print(fpath_out, 'saved.')
    
    def write_file(msg, output_path='./'):
        try:
            os.remove(output_path)
        except OSError:
            pass
    
        with open(output_path, 'w') as f:
            f.write(msg)
            print(output_path, 'saved.')
    
    
    def append_hgt_to_coords(coords, heights):
        coords2 = []
        for i in range(len(coords)):
            for hgt_m in heights:
                coords2.append(coords[i] + (hgt_m,))
    
        n_obs = len(coords2)
        print('effective number of observations (with vertical levels):', n_obs,
              ' on each level:', len(coords))
        return coords2
    
    
    def preamble(n_obs_3d_total, list_kinds):
        """Writes the header of an obs_seq.out file
        """
        lines_obstypedef = ''
        for kind in list_kinds:
            lines_obstypedef += '\n         '+str(obs_kind_nrs[kind])+' '+kind
    
        n_obs_str = str(n_obs_3d_total)
        num_obstypes = str(len(list_kinds))
        return """ obs_sequence
    obs_kind_definitions
               """+num_obstypes+' '+lines_obstypedef+"""
      num_copies:            0  num_qc:            0
      num_obs:            """+n_obs_str+"  max_num_obs:            "+n_obs_str+"""
      first:            1  last:            """+n_obs_str
    
    
    def determine_vert_coords(sat_channel, kind, obscfg):
        if not sat_channel:
            if 'SYNOP' in kind:
                vert_coord_sys = "-1"  # meters AGL
                vert_coords = [2, ]
            else:
                vert_coord_sys = "3"  # meters AMSL
                vert_coords = obscfg['heights']
        else:
            vert_coord_sys = "-2"  # undefined height
            vert_coords = ["-888888.0000000000", ]
        return vert_coord_sys, vert_coords
    
    
    def write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt):
        """Writes metadata str for an observation inside obs_seq.out
        """
        if sat_channel:
            sun_az = str(get_azimuth(lat0, lon0, time_dt))
            sun_zen = str(90. - get_altitude(lat0, lon0, time_dt))
            print('sunzen', sun_zen, 'sunazi', sun_az)
    
            return """visir
        """+sat_az+"""        """+sat_zen+"""        """+sun_az+"""
        """+sun_zen+"""
                12           4          21           """+str(sat_channel)+"""
        -888888.000000000
                1"""
        else:
            return ''
    
    
    def write_section(obs, last=False):
        """Returns the str of one observation inside an obs_seq.out file
    
        Args:
            obs (object)
            last (bool): True if this is the last observation in the obs_seq file
        """
        lon_rad = str(degr_to_rad(obs['lon']))
        lat_rad = str(degr_to_rad(obs['lat']))
    
        if last:
            line_link = "          "+str(obs['i']-1)+"           -1          -1"
        else:
            line_link = "        -1           "+str(obs['i']+1)+"          -1"
        
        return """
     OBS            """+str(obs['i'])+"""
    """+line_link+"""
    obdef
    loc3d
         """+lon_rad+"        "+lat_rad+"        "+str(obs['vert_coord'])+"     "+obs['vert_coord_sys']+"""
    kind
             """+str(obs['kind_nr'])+"""
    """+obs['appendix']+"""
    """+obs['secs_thatday']+"""     """+obs['dart_date_day']+"""
    """+str(obs['obserr_var'])
    
    
    def create_obs_seq_in(time_dt, list_obscfg, 
                                 output_path=cluster.dartrundir+'/obs_seq.in'):
        """Create obs_seq.in with multiple obs types in one file
    
        Args:
            time_dt (dt.datetime): time of observation
            list_obscfg (list of dict) :    configuration for observation types
                must have keys: 
                    - n_obs (int) : number of observations (must be a square of an integer: 4, 9, 1000, ...)
                    - obs_locations (str or tuple) in ['square_array_from_domaincenter', 'square_array_evenly_on_grid', ] 
                                                    or list of (lat, lon) coordinate tuples, in degrees north/east
                    - error_generate (float)
                    - error_assimilate (float or False) : False -> parameterized
                    - cov_loc_radius_km (float)
    
            obs_errors (np.array): contains observation errors, one for each observation
        """
        print('creating obs_seq.in:')
        time_dt = add_timezone_UTC(time_dt)
        dart_date_day, secs_thatday = get_dart_date(time_dt)
    
        txt = ''
        i_obs_total = 0
        for i_cfg, obscfg in enumerate(list_obscfg):
            n_obs = obscfg['n_obs']
    
            if obscfg['obs_locations'] == 'square_array_from_domaincenter':
                coords = col.square_array_from_domaincenter(n_obs, 
                                                            obscfg['distance_between_obs_km'])  # <---- must have variable
    
            elif obscfg['obs_locations'] == 'square_array_evenly_on_grid':
                coords = col.evenly_on_grid(n_obs)
    
            else:  # obs_locations given by iterable
                coords = obscfg['obs_locations']
    
            assert len(coords) == n_obs, (len(coords), n_obs)  # check if functions did what they supposed to
            for lat, lon in coords:
                assert (lat < 90) & (lat > -90)
                assert (lon < 180) & (lon > -180)
    
            kind = obscfg['kind']
            print('obstype', kind)
            sat_channel = obscfg.get('sat_channel', False)
    
            # add observation locations in the vertical at different levels
            vert_coord_sys, vert_coords = determine_vert_coords(sat_channel, kind, obscfg)
            coords = append_hgt_to_coords(coords, vert_coords)
            n_obs_3d_thistype = len(coords)
    
            # user defined generation error
            obserr_std = np.zeros(n_obs_3d_thistype) + obscfg["error_generate"]
    
            sat_info = write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt)
    
            for i_obs in range(n_obs_3d_thistype):
                i_obs_total += 1
                last = False
    
                is_last_obs_in_type = (i_obs == int(n_obs_3d_thistype)-1)
                is_last_obstype = (i_cfg == len(list_obscfg)-1)
    
                if is_last_obs_in_type and is_last_obstype:
                    last = True
    
                txt += write_section(dict(i=i_obs_total,
                                        kind_nr=obs_kind_nrs[kind],
                                        dart_date_day=dart_date_day,
                                        secs_thatday=secs_thatday,
                                        lon=coords[i_obs][1],
                                        lat=coords[i_obs][0],
                                        vert_coord=coords[i_obs][2],
                                        vert_coord_sys=vert_coord_sys,
                                        obserr_var=obserr_std[i_obs]**2,
                                        appendix=sat_info),
                                    last=last)
                
        n_obs_total = i_obs_total
        list_kinds = [a['kind'] for a in list_obscfg]
        pretxt = preamble(n_obs_total, list_kinds)
    
        alltxt = pretxt + txt
        write_file(alltxt, output_path=output_path)
    
    
    if __name__ == '__main__':
        # for testing
        time_dt = dt.datetime(2008, 7, 30, 13, 0)
    
    
        vis = dict(plotname='VIS 0.6µm',  plotunits='[1]',
                kind='MSG_4_SEVIRI_BDRF', sat_channel=1, 
                n_obs=961, obs_locations='square_array_evenly_on_grid',
                error_generate=0.03, error_assimilate=0.06,
                cov_loc_radius_km=32)
    
        # wv62 = dict(plotname='Brightness temperature WV 6.2µm', plotunits='[K]',
        #             kind='MSG_4_SEVIRI_TB', sat_channel=5, 
        #             n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
        #             error_generate=1., error_assimilate=False, 
        #             cov_loc_radius_km=20)
    
        # wv73 = dict(plotname='Brightness temperature WV 7.3µm',  plotunits='[K]',
        #             kind='MSG_4_SEVIRI_TB', sat_channel=6, 
        #             n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
        #             error_generate=1., error_assimilate=False,
        #             cov_loc_radius_km=20)
    
        # ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
        #             kind='MSG_4_SEVIRI_TB', sat_channel=9, 
        #             n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
        #             error_generate=5., error_assimilate=10.,
        #             cov_loc_radius_km=32)
    
        radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
                    kind='RADAR_REFLECTIVITY',             
                    n_obs=1, obs_locations=[(45,0),],
                    error_generate=2.5, error_assimilate=5.,
                    heights=np.arange(1000, 15001, 1000),
                    cov_loc_radius_km=20, cov_loc_vert_km=4)
    
        # t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
        #         kind='SYNOP_TEMPERATURE',             
        #         n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
        #         error_generate=0.1, error_assimilate=1.,
        #         cov_loc_radius_km=20, cov_loc_vert_km=3)
    
        # psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
        #             kind='SYNOP_SURFACE_PRESSURE', 
        #             n_obs=n_obs, obs_locations='square_array_evenly_on_grid',
        #             error_generate=50., error_assimilate=100.,
        #             cov_loc_radius_km=32, cov_loc_vert_km=5)
    
    
        create_obs_seq_in(time_dt, [radar])
    
        if False:
            error_assimilate = 5.*np.ones(n_obs*len(radar['heights']))
            import assim_synth_obs as aso
            aso.replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)