Working With Data From The CROCUS Doppler LIDAR Deployed During AGES+#
The CROCUS Scanning Doppler LIDAR was deployed to the NEIU Level Three rooftop node for July and August for NASA/NOAA AGES+. We have been running a crude ingest for the verticaly pointing data. This notebook gives an example of how to use and access this data.
For enquires contact Scott Collis scollis@anl.gov Paytsar Muradyan pmuradyan@anl.gov Max Grover mgrover@anl.gov
import numpy as np
from netCDF4 import Dataset
import os
import datetime
import xarray as xr
import pandas as pd
import matplotlib.dates as mdates
import matplotlib.dates as mdates
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
import pyart #This will be replaced by cmweather in the future
import act
import glob
%matplotlib inline
## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar data. Py-ART is partly
## supported by the U.S. Department of Energy as part of the Atmospheric
## Radiation Measurement (ARM) Climate Research Facility, an Office of
## Science user facility.
##
## If you use this software to prepare a publication, please cite:
##
## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[1], line 11
9 from matplotlib import pyplot as plt
10 from matplotlib.colors import LogNorm
---> 11 import pyart #This will be replaced by cmweather in the future
12 import act
13 import glob
File ~/miniconda3/envs/instrument-cookbooks-dev/lib/python3.10/site-packages/pyart/__init__.py:29
25 import importlib.metadata as _importlib_metadata
27 # import subpackages
28 # print out helpful message if build fails or importing from source tree
---> 29 from . import (
30 __check_build, # noqa
31 aux_io, # noqa
32 bridge, # noqa
33 config, # noqa
34 core, # noqa
35 correct, # noqa
36 filters, # noqa
37 graph, # noqa
38 io, # noqa
39 map, # noqa
40 retrieve, # noqa
41 testing, # noqa
42 util, # noqa
43 xradar, # noqa
44 )
45 from ._debug_info import _debug_info # noqa
47 # root level functions
File ~/miniconda3/envs/instrument-cookbooks-dev/lib/python3.10/site-packages/pyart/xradar/__init__.py:1
----> 1 from .accessor import Xradar, Xgrid # noqa
File ~/miniconda3/envs/instrument-cookbooks-dev/lib/python3.10/site-packages/pyart/xradar/accessor.py:8
1 """
2 Utilities for interfacing between xradar and Py-ART
3
4 """
6 import copy
----> 8 import datatree
9 import numpy as np
10 import pandas as pd
File ~/miniconda3/envs/instrument-cookbooks-dev/lib/python3.10/site-packages/datatree/__init__.py:2
1 # import public API
----> 2 from .datatree import DataTree
3 from .extensions import register_datatree_accessor
4 from .io import open_datatree
File ~/miniconda3/envs/instrument-cookbooks-dev/lib/python3.10/site-packages/datatree/datatree.py:33
31 from xarray.core.merge import dataset_update_method
32 from xarray.core.options import OPTIONS as XR_OPTS
---> 33 from xarray.core.utils import (
34 Default,
35 Frozen,
36 HybridMappingProxy,
37 _default,
38 either_dict_or_kwargs,
39 maybe_wrap_array,
40 )
41 from xarray.core.variable import Variable
43 from . import formatting, formatting_html
ImportError: cannot import name 'HybridMappingProxy' from 'xarray.core.utils' (/home/runner/miniconda3/envs/instrument-cookbooks-dev/lib/python3.10/site-packages/xarray/core/utils.py)
Open the dataset - This will be replaced with a THREDDS call
ds_ingested = act.io.read_arm_netcdf('/Volumes/T7/data/DL_ingested/crocus-neiu-dlvpt-a1-202308020000.nc')
Lets see what we have!
ds_ingested
/Users/scollis/miniconda3/envs/crocus/lib/python3.11/site-packages/xarray/core/formatting_html.py:22: DeprecationWarning: read_binary is deprecated. Use files() instead. Refer to https://importlib-resources.readthedocs.io/en/latest/using.html#migrating-from-legacy for migration advice.
read_binary(package, resource).decode("utf-8")
/Users/scollis/miniconda3/envs/crocus/lib/python3.11/site-packages/xarray/core/formatting_html.py:22: DeprecationWarning: read_binary is deprecated. Use files() instead. Refer to https://importlib-resources.readthedocs.io/en/latest/using.html#migrating-from-legacy for migration advice.
read_binary(package, resource).decode("utf-8")
<xarray.Dataset>
Dimensions: (range: 500, time: 144257)
Coordinates:
* range (range) float64 12.0 36.0 60.0 ... 1.196e+04 1.199e+04
* time (time) datetime64[ns] 2023-08-02T00:00:44 ... 2023-08-02...
Data variables:
radial_velocity (range, time) float64 dask.array<chunksize=(500, 144257), meta=np.ndarray>
beta (range, time) float64 dask.array<chunksize=(500, 144257), meta=np.ndarray>
intensity (range, time) float64 dask.array<chunksize=(500, 144257), meta=np.ndarray>
Attributes:
conventions: CF 1.10
site_ID: NEIU
CAMS_tag: CMS-SDL-001
datastream: neiu-dopplerlidar-vpt-a1
datalevel: a1
latitude: 41.9804526
longitude: -87.7196038
_file_dates: ['20230802']
_file_times: ['000044']
_datastream: neiu-dopplerlidar-vpt-a1
_arm_standards_flag: 1ds_ingested.intensity
<xarray.DataArray 'intensity' (range: 500, time: 144257)> dask.array<open_dataset-0c64766fbbc55c4b84caad558b84b479intensity, shape=(500, 144257), dtype=float64, chunksize=(500, 144257), chunktype=numpy.ndarray> Coordinates: * range (range) float64 12.0 36.0 60.0 ... 1.194e+04 1.196e+04 1.199e+04 * time (time) datetime64[ns] 2023-08-02T00:00:44 ... 2023-08-02T23:58:42
ds_ingested.intensity.plot()
<matplotlib.collections.QuadMesh at 0x1a32d0990>
Hmmm There is some data there but a a lot of data gathered around “6”
#time ranges
time1 = pd.Timestamp("2023-08-02 01:00:00")
time2 = pd.Timestamp("2023-08-02 02:00:00")
#Top of atmosphere
toa = 4000.0
#What ranges for Doppler Velocity
ds_ingested.intensity.sel(time=slice(time1,time2)).plot.pcolormesh(cmap=pyart.graph.cm_colorblind.ChaseSpectral)
my_axis = plt.gca()
my_axis.set_ylim([0,toa])
(0.0, 4000.0)
Lets look at possible values… Lets do a histogram
ds_ingested.intensity.plot.hist()
plt.yscale("log")
ok! So looking into the data for some reason the DL reports as Intensity + 1… likely to avoid issues with negative logs… so lets simply subtract one and scale the colormap to log
#time ranges
time1 = pd.Timestamp("2023-08-02 01:00:00")
time2 = pd.Timestamp("2023-08-02 02:00:00")
#Top of atmosphere
toa = 4000.0
#What ranges for Doppler Velocity
vrange = 5
fig, axs = plt.subplots( ncols=1, nrows=2, figsize=[20,8], constrained_layout=True)
#Units in the DL data are odd..
(ds_ingested.intensity -1.).sel(time=slice(time1,time2)).plot.pcolormesh(norm=LogNorm(vmax=3, vmin=0.001),
cmap=pyart.graph.cm_colorblind.ChaseSpectral, ax=axs[0])
plt.ylabel("Height")
axs[0].set_ylim([0,toa])
ds_ingested.radial_velocity.sel(time=slice(time1,time2)).plot.pcolormesh( cmap=pyart.graph.cm_colorblind.balance,
vmin=-vrange, vmax=vrange, ax=axs[1])
axs[1].set_ylim([0,toa])
(0.0, 4000.0)
Lovely… now how about some insight#
We know PBL heights are a key ask of the thematic groups. So lets play with looking at aggreating temporal data
This is inspired by the most excellent workshop run at the 2023 ARM-ASR PI meeting See: https://github.com/ARM-Development/ARM-Notebooks/blob/main/Tutorials/arm-asr-pi-meeting-2023/pandas_xarray_tutorial/Pandas_Xarray_intro.ipynb
Lets look at stats on an hour by hour window
window = 'time.hour'
#NB: For non standard intervals look at resample functions from the same ARM ASR Tutorial
new_ds = ds_ingested.groupby(window).mean()
(new_ds.intensity-1.).transpose().plot()
my_axis = plt.gca()
my_axis.set_ylim([0,toa])
(0.0, 4000.0)
new_ds.radial_velocity.transpose().plot(vmin=-2, vmax=2, cmap=pyart.graph.cm_colorblind.balance)
my_axis = plt.gca()
my_axis.set_ylim([0,toa])
(0.0, 4000.0)
Cool!#
But now lets look at standard deviation over the hour…
new_ds_stdev = ds_ingested.groupby(window).std()
new_ds_stdev.radial_velocity.transpose().plot()
my_axis = plt.gca()
my_axis.set_ylim([0,toa])
(0.0, 4000.0)
Ok… there is something there but the signal is swamped by the “texture” from when we are below SNR thresholds.
new_ds_stdev.radial_velocity.transpose().plot(vmax=3)
my_axis = plt.gca()
my_axis.set_ylim([0,toa])
(0.0, 4000.0)
The idea: Conditionally sample based on intensity…#
new_ds.intensity.plot.hist()
plt.yscale("log")
Woo Ho! Now lets make a data array which sets STDEV to zero when we think we are below detection limits…
filtered_da = new_ds_stdev.radial_velocity.where(new_ds.intensity > 1.01, 0)
filtered_da.transpose().plot(vmax=2)
my_axis = plt.gca()
my_axis.set_ylim([0,toa])
(0.0, 4000.0)
#This will take a LONG time to run
#Top of atmosphere
toa = 4000.0
#What ranges for Doppler Velocity
vrange = 5
fig, axs = plt.subplots( ncols=1, nrows=3, figsize=[20,12], constrained_layout=True)
#Units in the DL data are odd..
(ds_ingested.intensity -1.).plot.pcolormesh(norm=LogNorm(vmax=3, vmin=0.001),
cmap=pyart.graph.cm_colorblind.ChaseSpectral, ax=axs[0])
plt.ylabel("Height")
axs[0].set_ylim([0,toa])
ds_ingested.radial_velocity.plot.pcolormesh( cmap=pyart.graph.cm_colorblind.balance,
vmin=-vrange, vmax=vrange, ax=axs[1])
axs[1].set_ylim([0,toa])
filtered_da.transpose().plot.pcolormesh(vmax=2, ax=axs[2])
axs[2].set_ylim([0,toa])