Skip to content
Snippets Groups Projects
Select Git revision
  • 66c0149f62aeced19f38eaa228c003903f28162c
  • consistent_config default protected
2 results

test.ipynb

Blame
  • calculate_obs_locations.py 3.73 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
    import xarray as xr
    
    from exp_config import exp
    from server_config import cluster
    
    #####################
    # Global variables
    
    # position on Earth for DART, domain center when coords_from_domaincenter=True
    lat0_center = 45.
    lon0_center = 0.
    
    radius_earth_meters = 6.371*1E6
    
    def square_array_from_domaincenter(n_obs, distance_between_obs_km):
        """
        Create equally spaced grid for satellite observations every 4 km
        e.g. ny,nx = 10
        -> obs locations are from -5 to +5 times dy in south_north direction
        and from -5 to +5 times dx in west_east direction
    
        Returns
            tuple of (lat, lon) coordinates
        """
        coords = []
        nx, ny = int(np.sqrt(n_obs)), int(np.sqrt(n_obs))
    
        m_per_degree = 2*np.pi*radius_earth_meters/360
        distance_between_obs_meters = distance_between_obs_km*1000.
        dy_4km_in_degree = distance_between_obs_meters/m_per_degree
    
        for iy in range(int(-ny/2), int(ny/2+1)):
            for ix in range(int(-nx/2), int(nx/2+1)):
    
                lat = lat0_center + iy*dy_4km_in_degree
                m_per_degree_x = 2*np.pi*radius_earth_meters*np.sin(lat/180*np.pi)/360
                dx_4km_in_degree = distance_between_obs_meters/m_per_degree_x
                lon = lon0_center + ix*dx_4km_in_degree
                coords.append((lat, lon))
        return coords
    
    def evenly_on_grid(n_obs, omit_covloc_radius_on_boundary=True):
        """Observations spread evenly over domain
        
        omit_covloc_radius_on_boundary : leave out a distance to the border of the domain
                                         so that the assimilation increments are zero on the boundary
                                         distance to boundary = 50 km
    
        Returns
            tuple of (lat, lon) coordinates
        """
        fcoords = cluster.geo_em_for_WRF_ideal  # TODO: in case of WRF real, need to find a file providing coordinates, e.g. nature wrfout
        ds = xr.open_dataset(fcoords)
    
        lons = ds.XLONG_M.isel(Time=0).values
        lats = ds.XLAT_M.isel(Time=0).values
        n_obs_x = int(np.sqrt(n_obs))
        nx = len(ds.south_north)  # number of gridpoints in one direction
        model_dx_km = exp.model_dx/1000
        print('assuming', model_dx_km, 'km model grid')
    
    
        if omit_covloc_radius_on_boundary:  #  in order to avoid an increment step on the boundary
            skip_km = 50  
            skip_gridpoints = int(skip_km/model_dx_km)  # skip this many gridpoints on each side
    
            gridpoints_left = nx - 2*skip_gridpoints
            # now spread observations evenly across the space left
            gridpoints_between_obs = int(gridpoints_left/(n_obs_x-1))
        else:
            gridpoints_between_obs = int(nx/n_obs_x)  # number of gridpoints / number of observations in one direction
            skip_gridpoints = int(gridpoints_between_obs/2)    # to center the observations in the domain
        
        km_between_obs = gridpoints_between_obs*model_dx_km
        print('observation density: one observation every', km_between_obs, 'km /', 
        gridpoints_between_obs,'gridpoints \n', 'leaving a domain boundary on each side of', 
        skip_gridpoints, 'gridpoints or', skip_gridpoints*model_dx_km, 'km')
        # skip_gridpoints = space we have to leave out on the boundary of the domain
        # in order to have the observations centered in the domain
    
        coords = []
        for i in range(n_obs_x):
            for j in range(n_obs_x):
                coords.append((lats[skip_gridpoints+i*gridpoints_between_obs, skip_gridpoints+j*gridpoints_between_obs],
                                lons[skip_gridpoints+i*gridpoints_between_obs, skip_gridpoints+j*gridpoints_between_obs]))
        return coords