Any ideas where I am going wrong merging two local radar files

[ahuang11]

I am trying to do something similar to:

!wget -nc https://nomads.ncep.noaa.gov/pub/data/nccf/radar/nexrad_level2/KILX/KILX_20210115_103151.bz2
!wget -nc https://nomads.ncep.noaa.gov/pub/data/nccf/radar/nexrad_level2/KLOT/KLOT_20210115_103558.bz2
import pyart
radar1 = pyart.io.read_nexrad_archive('KILX_20210115_103151.bz2')
radar2 = pyart.io.read_nexrad_archive('KLOT_20210115_103558.bz2')
grid = pyart.map.grid_from_radars(
    (radar1, radar2),
    grid_shape=(1, 201, 201),
    grid_limits=((1000, 1000), (-100000, 100000), (-100000, 100000)),
    grid_origin = (40, -88)
)
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(grid.fields['reflectivity']['data'][0],
          origin='lower', vmin=0, vmax=48)
plt.show()

But with OpenMosaic, so I cut and paste some stuff from (I edited stuff in if __name == "main"):
https://github.com/jthielen/OpenMosaic/blob/main/examples/process_full_conus.py

But I get all NaNs at the end.

import pyproj
import shapely
import pandas as pd
import numpy as np
import xarray as xr
from openmosaic.gridding import Gridder, generate_rectangular_grid
from openmosaic.utils import get_nexrad_sites
from dask.distributed import Client, LocalCluster
cluster = LocalCluster()
cluster.scale(6)
client = Client(cluster)

# Function definitions
def aggregate_to_dataset(fields, grid_params, cf_attrs, analysis_time, *subgrid_data_and_weights):
    """Aggregate results from individual radar subgrids onto main grid."""
    grid_shape = (grid_params['nz'], grid_params['ny'], grid_params['nx'], len(fields))
    grid_sum = np.zeros(grid_shape, dtype='float32')
    grid_wsum = np.zeros(grid_shape, dtype='float32')
    field_metadata = {field: '' for field in fields}  # If no valid data, have to mock metadata values
    for subgrid in subgrid_data_and_weights:
        if subgrid is not None:
            field_metadata = subgrid['field_metadata']
            grid_sum[:, subgrid['yi_min']:subgrid['yi_max'], subgrid['xi_min']:subgrid['xi_max'], :] += subgrid['sum']
            grid_wsum[:, subgrid['yi_min']:subgrid['yi_max'], subgrid['xi_min']:subgrid['xi_max'], :] += subgrid['wsum']

    # Apply the weighting and masking
    mweight = np.ma.masked_equal(grid_wsum, 0)
    msum = np.ma.masked_array(grid_sum, mweight.mask)
    grids = {f: (msum[..., i] / mweight[..., i]) for i, f in enumerate(fields)}

    # Create Dataset from dictionary
    ds = generate_rectangular_grid(
        grid_params['nx'],
        grid_params['ny'],
        grid_params['dx'],
        grid_params['dy'],
        cf_attrs,
        x0=grid_params['x_min'],
        y0=grid_params['y_min']
    )
    ds = ds.assign_coords(
        z=xr.DataArray(
            np.arange(grid_params['nz']) * grid_params['dz'] + grid_params['z_min'],
            dims=('z',),
            name='z',
            attrs={
                'standard_name': 'altitude',
                'units': 'meter',
                'positive': 'up'
            }
        ),
        time=xr.DataArray(
            analysis_time,
            name='time',
            attrs={
                'long_name': 'Time of radar analysis'
            }
        )
    )
    ds = ds.assign(
        {
            field: xr.Variable(
                ('z', 'y', 'x'),
                grids[field],
                field_metadata[field]
            )
            for field in grids
        }
    )

    return ds


def post_ops(dataset):
    """Post subbatch dataset operations to go from 3D to 2D.
    
    Notes
    -----
    Currently only supports column-maximum reflectivity calculations!
    """
    return dataset.max('z').drop_vars(['lon', 'lat'])



if __name__ == '__main__':
    !wget -nc https://nomads.ncep.noaa.gov/pub/data/nccf/radar/nexrad_level2/KILX/KILX_20210115_103151.bz2
    !wget -nc https://nomads.ncep.noaa.gov/pub/data/nccf/radar/nexrad_level2/KLOT/KLOT_20210115_103558.bz2
    !wget -nc https://opendata.arcgis.com/datasets/85c59fe951504e9b9919e24d7a684084_3.geojson

    nexrad_files = [
        'KILX_20210115_103151.bz2'
        'KLOT_20210115_103558.bz2'
    ]
    temp_reg_dir = '.'
    nexrad_file = '85c59fe951504e9b9919e24d7a684084_3.geojson'
    nexrad_site_ids = ['KILX', 'KLOT']
    timestamp = pd.Timestamp('20210115 10:35')
    sweep_interval = 228

    cf_attrs = {
        'grid_mapping_name': 'lambert_conformal_conic',
        'latitude_of_projection_origin': 40,
        'longitude_of_central_meridian': 272.5,
        'standard_parallel': 38.5,
        'earth_radius': 6371229.0
    }
    crs = pyproj.CRS.from_cf(cf_attrs)
    proj = pyproj.Proj(crs)
    radar_sites = get_nexrad_sites(nexrad_file, crs)
    radar_sites['x'] = radar_sites['proj_geom'].x
    radar_sites['y'] = radar_sites['proj_geom'].y
    radar_sites_no_geo = pd.DataFrame(radar_sites[['x', 'y']])

    g = Gridder(cf_attrs, nz=24, dz=1000., z_min=1000.)
    horizontal_params = {
        'x_min': -703215.388043841,
        'x_max': 703215.388043841,
        'y_min': -250000.1279402958,
        'y_max': 250000.1279402958
    }
    r_max = 3.0e7
    # Subset the radar definition file to only include those sites in your domain
    radar_search_area = shapely.geometry.box(
        horizontal_params['x_min'] - r_max,
        horizontal_params['y_min'] - r_max,
        horizontal_params['x_max'] + r_max,
        horizontal_params['y_max'] + r_max
    )
    g.assign_from_subbatch_and_spacing(horizontal_params, 2000.)
    g.cache_dir = temp_reg_dir
    g.prepare_radar_subgrids(set(nexrad_site_ids), radar_sites_no_geo, r_max)

    # Go from nexrad file on disk to 4D subgrid at analysis time (extra dim is fields)
    fields = ['reflectivity']
    subgrid_data_and_weights_futures = client.map(
        g.process_nexrad_file_to_subgrids_and_weights,
        nexrad_files,
        nexrad_site_ids,
        fields=fields,
        analysis_time=timestamp,
        sweep_interval=sweep_interval,
        priority=2
    )

    # Go from collection of subgrid data and weights to full 3D xarray dataset
    # def aggregate_to_dataset(*subgrid_data_and_weights)
    dataset_full = client.submit(
        aggregate_to_dataset,
        fields,
        g.grid_params,
        g.cf_attrs,
        timestamp,
        *subgrid_data_and_weights_futures,
        priority=3
    )

    datasets = [dataset_full]

    gathered_datasets = client.gather(datasets)
    output = xr.concat(gathered_datasets, dim='time')

    print(output['reflectivity'].max())

[jthielen]

Unfortunately, I can no longer access those radar files in order to reproduce your example, as they return 404 Not Found.

However, I suspect what is going on is that the warning messages in g.process_nexrad_file_to_subgrids_and_weights are being swallowed/ignored somewhere, as an all nan grid is an expected result if no valid radar data (either all missing files, all outside of domain, or all outside of sweep interval) are found. I'll put that on my to-do list to troubleshoot in the gridder class.

In case you were to go back and troubleshoot this on your own in the meantime with updated radar files, there are two bugs I'd recommend fixing in your modified example. First, you appear to be missing a comma in the nexrad_files list...instead of

nexrad_files = [
    'KILX_20210115_103151.bz2'
    'KLOT_20210115_103558.bz2'
]

(which creates a single item list with an invalid file name) it should be

nexrad_files = [
    'KILX_20210115_103151.bz2',
    'KLOT_20210115_103558.bz2'
]

Second, sweep_interval needs to be a time delta compatible with the pd.Timestamp, since time arithmetic is performed. I'd recommend a np.timedelta64 as the easiest choice, so sweep_interval = np.timedelta64(228, 's').

Hopefully that helps, and sorry for this not being that easy to use for your use case at the current moment!

[ahuang11]

Okay I got something to show...

  
import pyproj
import shapely
import pandas as pd
import numpy as np
import xarray as xr
from openmosaic.gridding import Gridder, generate_rectangular_grid
from openmosaic.utils import get_nexrad_sites
from dask.distributed import Client, LocalCluster
client = Client()

!wget -nc https://nomads.ncep.noaa.gov/pub/data/nccf/radar/nexrad_level2/KILX/KILX_20210119_080351.bz2
!wget -nc https://nomads.ncep.noaa.gov/pub/data/nccf/radar/nexrad_level2/KLOT/KLOT_20210119_075832.bz2

    
NEXRAD_FILES = [
    'KILX_20210119_080351.bz2',
    'KLOT_20210119_075832.bz2'
]

TIMESTAMP = pd.Timestamp('20210119 08:00')

# Function definitions
def aggregate_to_dataset(fields, grid_params, cf_attrs, analysis_time, *subgrid_data_and_weights):
    """Aggregate results from individual radar subgrids onto main grid."""
    grid_shape = (grid_params['nz'], grid_params['ny'], grid_params['nx'], len(fields))
    grid_sum = np.zeros(grid_shape, dtype='float32')
    grid_wsum = np.zeros(grid_shape, dtype='float32')
    field_metadata = {field: '' for field in fields}  # If no valid data, have to mock metadata values
    for subgrid in subgrid_data_and_weights:
        if subgrid is not None:
            field_metadata = subgrid['field_metadata']
            grid_sum[:, subgrid['yi_min']:subgrid['yi_max'], subgrid['xi_min']:subgrid['xi_max'], :] += subgrid['sum']
            grid_wsum[:, subgrid['yi_min']:subgrid['yi_max'], subgrid['xi_min']:subgrid['xi_max'], :] += subgrid['wsum']

    # Apply the weighting and masking
    mweight = np.ma.masked_equal(grid_wsum, 0)
    msum = np.ma.masked_array(grid_sum, mweight.mask)
    grids = {f: (msum[..., i] / mweight[..., i]) for i, f in enumerate(fields)}

    # Create Dataset from dictionary
    ds = generate_rectangular_grid(
        grid_params['nx'],
        grid_params['ny'],
        grid_params['dx'],
        grid_params['dy'],
        cf_attrs,
        x0=grid_params['x_min'],
        y0=grid_params['y_min']
    )
    ds = ds.assign_coords(
        z=xr.DataArray(
            np.arange(grid_params['nz']) * grid_params['dz'] + grid_params['z_min'],
            dims=('z',),
            name='z',
            attrs={
                'standard_name': 'altitude',
                'units': 'meter',
                'positive': 'up'
            }
        ),
        time=xr.DataArray(
            analysis_time,
            name='time',
            attrs={
                'long_name': 'Time of radar analysis'
            }
        )
    )
    ds = ds.assign(
        {
            field: xr.Variable(
                ('z', 'y', 'x'),
                grids[field],
                field_metadata[field]
            )
            for field in grids
        }
    )

    return ds


def post_ops(dataset):
    """Post subbatch dataset operations to go from 3D to 2D.
    
    Notes
    -----
    Currently only supports column-maximum reflectivity calculations!
    """
    return dataset.max('z').drop_vars(['lon', 'lat'])



if __name__ == '__main__':
    !wget -nc https://opendata.arcgis.com/datasets/85c59fe951504e9b9919e24d7a684084_3.geojson

    temp_reg_dir = '.'
    nexrad_file = '85c59fe951504e9b9919e24d7a684084_3.geojson'
    nexrad_site_ids = ['KILX', 'KLOT']
    sweep_interval = np.timedelta64(228, 's')

    cf_attrs = {
        'grid_mapping_name': 'lambert_conformal_conic',
        'latitude_of_projection_origin': 40,
        'longitude_of_central_meridian': 272.5,
        'standard_parallel': 38.5,
        'earth_radius': 6371229.0
    }
    crs = pyproj.CRS.from_cf(cf_attrs)
    proj = pyproj.Proj(crs)
    radar_sites = get_nexrad_sites(nexrad_file, crs)
    radar_sites['x'] = radar_sites['proj_geom'].x
    radar_sites['y'] = radar_sites['proj_geom'].y
    radar_sites_no_geo = pd.DataFrame(radar_sites[['x', 'y']])

    g = Gridder(cf_attrs, nz=24, dz=1000., z_min=1000.)
    horizontal_params = {
        'x_min': -203215.388043841,
        'x_max': 50215.388043841,
        'y_min': 180000.1279402958,
        'y_max': 350000.1279402958
    }
    r_max = 3.0e7
    # Subset the radar definition file to only include those sites in your domain
    radar_search_area = shapely.geometry.box(
        horizontal_params['x_min'] - r_max,
        horizontal_params['y_min'] - r_max,
        horizontal_params['x_max'] + r_max,
        horizontal_params['y_max'] + r_max
    )
    g.assign_from_subbatch_and_spacing(horizontal_params, 2000.)
    g.cache_dir = temp_reg_dir
    g.prepare_radar_subgrids(set(nexrad_site_ids), radar_sites_no_geo, r_max)

    # Go from nexrad file on disk to 4D subgrid at analysis time (extra dim is fields)
    fields = ['reflectivity']
    subgrid_data_and_weights_futures = client.map(
        g.process_nexrad_file_to_subgrids_and_weights,
        NEXRAD_FILES,
        nexrad_site_ids,
        fields=fields,
        analysis_time=TIMESTAMP,
        sweep_interval=sweep_interval,
        priority=2
    )

    # Go from collection of subgrid data and weights to full 3D xarray dataset
    # def aggregate_to_dataset(*subgrid_data_and_weights)
    dataset_full = client.submit(
        aggregate_to_dataset,
        fields,
        g.grid_params,
        g.cf_attrs,
        TIMESTAMP,
        *subgrid_data_and_weights_futures,
        priority=3
    )

    datasets = [dataset_full]

    gathered_datasets = client.gather(datasets)
    output = xr.concat(gathered_datasets, dim='time')

    output['reflectivity'].isel(z=0).plot()

Is there a way to automatically figure out horizontal params? Or set it in degrees?

[jthielen]

Nice, glad to hear you got it figured out! And thanks for the example; it will be good to use as a guide towards where simplifications to the Gridder class would be helpful.

To automatically figure out horizontal params, you could use a pyproj transformer to convert between any coordinate systems that you'd want. Otherwise, you should also be able to easily grid onto a lat/lon grid instead of LCC. Just specify the cf_attrs to be

cf_attrs = {
    "grid_mapping_name": "latitude_longitude"
}

(with any optional globe param controls that you'd want), and then you just set the horizontal params in degrees (e.g., make sure dx and dy get changed to something like 0.02).