openradar / xradar

A tool to work in weather radar data in xarray
https://docs.openradarscience.org/projects/xradar
MIT License
98 stars 17 forks source link

Add an Example Using `open_mfdataset` with an angle preprocessing function #64

Closed mgrover1 closed 1 year ago

mgrover1 commented 1 year ago

Description

At the Pangeo working meeting on Friday, an xradar user suggested adding an example reading in multiple SIGMET files, applying the angle re-indexing, and reading into a single xarray object for the lowest sweep in the volume scan.

What I Did

It would look something like this:

import glob
import xradar

# Reindex function
def reindex_angles():
   ...

# List of files
sigmet_files = glob.glob('...')

# Read the data using preprocessing function
ds = xr.open_mfdataset(sigmet_files, engine='sigmet', preprocess=reindex_angles, concat_dim='volume')
kmuehlbauer commented 1 year ago

Just to understand that correctly, the layout of the resulting dataset should be:

Dimensions: (volume: 10, azimuth: 360, range: 1000)

for a file with 10 elevations, 360 azimuth rays and 1000 range bins?

Also for xr.open_mfdataset kwarg group needs to be given which refers to the wanted sweep. I think this doesn't work with a one line using xr.open_mfdataset.

It would be great if you could elaborate a bit more on the file layout and the final dataset layout to get the full picture.

mgrover1 commented 1 year ago

Here is a better illustration of the code! This worked... with each file including a single sweep.


import glob
import xradar

def fix_angle(ds):
    ds['time'] = ds.time.load() # Convert time from dask to numpy
    angle_dict = xd.util.extract_angle_parameters(ds)
    start_ang = 0 # Set consistent start/end values
    stop_ang = 360
    angle_res = 1.
    direction = angle_dict["direction"]

    # first find exact duplicates and remove
    ds = xd.util.remove_duplicate_rays(ds)

    # second reindex according to retrieved parameters
    ds = xd.util.reindex_angle(
        ds, start_ang, stop_ang, angle_res, direction, method="nearest"
    )

    ds = ds.expand_dims('volume_time') # Expand for volumes for concatenation

    ds['volume_time'] = [ds.time.min('azimuth').values]

    return ds

# List of files
sigmet_files = glob.glob('...')

# Concatenate in xarray ds
ds = xr.open_mfdataset(fns,
                       preprocess=fix_angle,
                       engine='iris',
                       group='sweep_0',
                       concat_dim='volume_time',
                       combine='nested')
kmuehlbauer commented 1 year ago

@mgrover1 Looks great so far. This could be expanded to read timeseries of volumes, too.

mgrover1 commented 1 year ago

Yup! Thanks! I'll draft up a PR this week.