Skip to article frontmatterSkip to article content

Simulation vs Observational Data of Shallow Cumulus Clouds over the Southern Great Plains on April 4th, 2019

Import

# Lasso Simulation Data

# import dask
from datetime import datetime
from distributed import Client
import numpy as np
import xarray as xr
import xwrf
import s3fs
import fsspec
import xarray as xr
import glob

import matplotlib.pyplot as plt

Spin up a Dask Cluster

We will use Dask here to access the data in a parallel/distributed manner.

client = Client()
client
Loading...

Access LASSO SGP Data from the NSF Jetstream Cloud

A subset of the LASSO Shallow Cumulus Experiment over the Southern Great Plains site has been made available on a cloud bucket, hosted through Project Pythia. These datasets were originally accessed through the LASSO bundle browser, untarred, then uploaded to the cloud bucket. We focus exclusively on the April 4, 2019 case.

Data were obtained from the Atmospheric Radiation Measurement (ARM) Program sponsored by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, Climate and Environmental Sciences Division.

Set the URL to access the data

The data are stored on a bucket, which is a web-accessible place where we can remotely stream the data, without downloading directly. The bucket is located on the NSF jetstream cloud, which we can see below. We then use fsspec to easily list the directrories and load in the data.

Below we set the url, then list (glob) the directories in the bucket.

# Set the URL and path for the cloud
URL = 'https://js2.jetstream-cloud.org:8001/'
path = f'pythia/lasso-sgp'


fs = fsspec.filesystem("s3", anon=True, client_kwargs=dict(endpoint_url=URL))

fs.glob(f"{path}/*")
['pythia/lasso-sgp/high_res_obs', 'pythia/lasso-sgp/metrics', 'pythia/lasso-sgp/sim0001', 'pythia/lasso-sgp/sim0002', 'pythia/lasso-sgp/sim0003', 'pythia/lasso-sgp/sim0004', 'pythia/lasso-sgp/sim0005', 'pythia/lasso-sgp/sim0006', 'pythia/lasso-sgp/sim0007', 'pythia/lasso-sgp/sim0008']

We notice that there are 8 simulations, as well as observations in the bucket. We are going to start with the fourth simulation, setting a path to the actual output, which is under /raw_model/

case_date = datetime(2019, 4, 4)
sim_id = 4

# Read the wrfstat files
wrfstat_pattern = f's3://{path}/sim000{sim_id}/raw_model/wrfstat*'

# Read the wrfout files
wrfout_pattern = f's3://{path}/sim000{sim_id}/raw_model/wrfout*'

wrfstat_files = sorted(fs.glob(wrfstat_pattern))
wrfout_files = sorted(fs.glob(wrfout_pattern))

Now that we have lists of files, we setup a path to read into xarray since we need the bucket information as well.

wrfstat_file_list = [fs.open(file) for file in wrfstat_files]
wrfout_file_list = [fs.open(file) for file in wrfout_files]

Load Data Using Xarray and View Variables

We have a single WRF stat file, which we can load into xarray, then postprocess with xwrf.

ds_stat = xr.open_mfdataset(wrfstat_file_list, engine='h5netcdf').xwrf.postprocess()
ds_stat
---------------------------------------------------------------------------
ClientError                               Traceback (most recent call last)
File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:114, in _error_wrapper(func, args, kwargs, retries)
    113 try:
--> 114     return await func(*args, **kwargs)
    115 except S3_RETRYABLE_ERRORS as e:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/aiobotocore/context.py:36, in with_current_context.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
     35     await resolve_awaitable(hook())
---> 36 return await func(*args, **kwargs)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/aiobotocore/client.py:424, in AioBaseClient._make_api_call(self, operation_name, api_params)
    423     error_class = self.exceptions.from_code(error_code)
--> 424     raise error_class(parsed_response, operation_name)
    425 else:

ClientError: An error occurred (PreconditionFailed) when calling the GetObject operation: None

The above exception was the direct cause of the following exception:

OSError                                   Traceback (most recent call last)
File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:2390, in S3File._fetch_range(self, start, end)
   2389 try:
-> 2390     return _fetch_range(
   2391         self.fs,
   2392         self.bucket,
   2393         self.key,
   2394         self.version_id,
   2395         start,
   2396         end,
   2397         req_kw=self.req_kw,
   2398     )
   2400 except OSError as ex:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:2559, in _fetch_range(fs, bucket, key, version_id, start, end, req_kw)
   2558 logger.debug("Fetch: %s/%s, %s-%s", bucket, key, start, end)
-> 2559 return sync(fs.loop, _inner_fetch, fs, bucket, key, version_id, start, end, req_kw)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/fsspec/asyn.py:103, in sync(loop, func, timeout, *args, **kwargs)
    102 elif isinstance(return_result, BaseException):
--> 103     raise return_result
    104 else:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/fsspec/asyn.py:56, in _runner(event, coro, result, timeout)
     55 try:
---> 56     result[0] = await coro
     57 except Exception as ex:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:2577, in _inner_fetch(fs, bucket, key, version_id, start, end, req_kw)
   2575         resp["Body"].close()
-> 2577 return await _error_wrapper(_call_and_read, retries=fs.retries)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:146, in _error_wrapper(func, args, kwargs, retries)
    145 err = translate_boto_error(err)
--> 146 raise err

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:114, in _error_wrapper(func, args, kwargs, retries)
    113 try:
--> 114     return await func(*args, **kwargs)
    115 except S3_RETRYABLE_ERRORS as e:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:2564, in _inner_fetch.<locals>._call_and_read()
   2563 async def _call_and_read():
-> 2564     resp = await fs._call_s3(
   2565         "get_object",
   2566         Bucket=bucket,
   2567         Key=key,
   2568         Range="bytes=%i-%i" % (start, end - 1),
   2569         **version_id_kw(version_id),
   2570         **req_kw,
   2571     )
   2572     try:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:371, in S3FileSystem._call_s3(self, method, *akwarglist, **kwargs)
    370 additional_kwargs = self._get_s3_method_kwargs(method, *akwarglist, **kwargs)
--> 371 return await _error_wrapper(
    372     method, kwargs=additional_kwargs, retries=self.retries
    373 )

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:146, in _error_wrapper(func, args, kwargs, retries)
    145 err = translate_boto_error(err)
--> 146 raise err

OSError: [Errno 22] None

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 ds_stat = xr.open_mfdataset(wrfstat_file_list, engine='h5netcdf').xwrf.postprocess()
      2 ds_stat

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/xarray/backends/api.py:1634, in open_mfdataset(paths, chunks, concat_dim, compat, preprocess, engine, data_vars, coords, combine, parallel, join, attrs_file, combine_attrs, errors, **kwargs)
   1632 for p in paths1d:
   1633     try:
-> 1634         ds = open_(p, **open_kwargs)
   1635         datasets.append(ds)
   1636     except Exception as e:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/xarray/backends/api.py:606, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, create_default_indexes, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
    594 decoders = _resolve_decoders_kwargs(
    595     decode_cf,
    596     open_backend_dataset_parameters=backend.open_dataset_parameters,
   (...)    602     decode_coords=decode_coords,
    603 )
    605 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 606 backend_ds = backend.open_dataset(
    607     filename_or_obj,
    608     drop_variables=drop_variables,
    609     **decoders,
    610     **kwargs,
    611 )
    612 ds = _dataset_from_backend_dataset(
    613     backend_ds,
    614     filename_or_obj,
   (...)    625     **kwargs,
    626 )
    627 return ds

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/xarray/backends/h5netcdf_.py:540, in H5netcdfBackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, format, group, lock, invalid_netcdf, phony_dims, decode_vlen_strings, driver, driver_kwds, storage_options)
    537 emit_phony_dims_warning, phony_dims = _check_phony_dims(phony_dims)
    539 filename_or_obj = _normalize_filename_or_obj(filename_or_obj)
--> 540 store = H5NetCDFStore.open(
    541     filename_or_obj,
    542     format=format,
    543     group=group,
    544     lock=lock,
    545     invalid_netcdf=invalid_netcdf,
    546     phony_dims=phony_dims,
    547     decode_vlen_strings=decode_vlen_strings,
    548     driver=driver,
    549     driver_kwds=driver_kwds,
    550     storage_options=storage_options,
    551 )
    553 store_entrypoint = StoreBackendEntrypoint()
    555 ds = store_entrypoint.open_dataset(
    556     store,
    557     mask_and_scale=mask_and_scale,
   (...)    563     decode_timedelta=decode_timedelta,
    564 )

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/xarray/backends/h5netcdf_.py:201, in H5NetCDFStore.open(cls, filename, mode, format, group, lock, autoclose, invalid_netcdf, phony_dims, decode_vlen_strings, driver, driver_kwds, storage_options)
    198     source.getvalue = filename.getbuffer
    200 if isinstance(filename, io.IOBase) and mode == "r":
--> 201     magic_number = read_magic_number_from_file(filename)
    202     if not magic_number.startswith(b"\211HDF\r\n\032\n"):
    203         raise ValueError(
    204             f"{magic_number!r} is not the signature of a valid netCDF4 file"
    205         )

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/xarray/core/utils.py:775, in read_magic_number_from_file(filename_or_obj, count)
    773 if filename_or_obj.tell() != 0:
    774     filename_or_obj.seek(0)
--> 775 magic_number = filename_or_obj.read(count)
    776 filename_or_obj.seek(0)
    777 return magic_number

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/fsspec/spec.py:2122, in AbstractBufferedFile.read(self, length)
   2119 if length == 0:
   2120     # don't even bother calling fetch
   2121     return b""
-> 2122 out = self.cache._fetch(self.loc, self.loc + length)
   2124 logger.debug(
   2125     "%s read: %i - %i %s",
   2126     self,
   (...)   2129     self.cache._log_stats(),
   2130 )
   2131 self.loc += len(out)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/fsspec/caching.py:287, in ReadAheadCache._fetch(self, start, end)
    285 end = min(self.size, end + self.blocksize)
    286 self.total_requested_bytes += end - start
--> 287 self.cache = self.fetcher(start, end)  # new block replaces old
    288 self.start = start
    289 self.end = self.start + len(self.cache)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.14/site-packages/s3fs/core.py:2401, in S3File._fetch_range(self, start, end)
   2390     return _fetch_range(
   2391         self.fs,
   2392         self.bucket,
   (...)   2397         req_kw=self.req_kw,
   2398     )
   2400 except OSError as ex:
-> 2401     if ex.args[0] == errno.EINVAL and "pre-conditions" in ex.args[1]:
   2402         raise FileExpired(
   2403             filename=self.details["name"], e_tag=self.details.get("ETag")
   2404         ) from ex
   2405     else:

TypeError: argument of type 'NoneType' is not a container or iterable
# Plotting wrfstat variables...

# path_shcu_root = "/gpfs/wolf2/arm/atm124/world-shared/arm-summer-school-2024/lasso_tutorial/ShCu/untar/"  # on cumulus
path_shcu_root = "/data/project/ARM_Summer_School_2024_Data/lasso_tutorial/ShCu/untar"  # on Jupyter

case_date = datetime(2019, 4, 4)
sim_id = 4

ds_stat = xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id:04d}/raw_model/wrfstat_d01_{case_date:%Y-%m-%d_12:00:00}.nc")
ds_stat

Plot Variables and Modify as Desired

xwrf automatically corrected the time for us! So now we can focus on subsetting given an hour. In this case, we are interested in 1700 UTC.

ds_stat.CSV_LWC
hour_to_plot = 17

# Time series:
ds_stat["CST_LWP"].plot()
plt.show()

# Profile at a selected time (plots sideways, though, since we are being lazy):
ds_stat["CSP_LWC"].sel(Time=f"{case_date:%Y-%m-%d} {hour_to_plot}:00").plot()
plt.show()

# X-Y slice for a selected time:
ds_stat["CSS_LWP"].sel(Time=f"{case_date:%Y-%m-%d} {hour_to_plot}:00").plot()
plt.show()

# A vertical slice from the volume at a selected time:
# We'll assign the vertical coordinate values for this one and hide the cloud-free upper atmosphere.
plot_data = ds_stat["CSV_LWC"].assign_coords(height=(ds_stat["CSP_Z"]))
plot_data.sel(Time=f"{case_date:%Y-%m-%d} {hour_to_plot}:00", y=-12450).plot(y="height", ylim=[0, 1500])
plt.show()

# Add lines and modify variables to plot desired figures... In this notebook, we plotted

Plot wind speed

Let’s start with the first 10 files from the WRF simulation.

# Note the extra details required by open_mfdataset to connect the files together in time.
ds_wrf = xr.open_mfdataset(wrfout_file_list[:10],
                           combine="nested",
                           concat_dim="Time")

Since we did not use xwrf, this time, let’s fix the times.

ds_wrf["Time"] = ds_wrf["XTIME"]

Visualize the Wind Vectors and Expose Destaggering

Next, we plot wind vectors at a selected level to demonstrate how to destagger the wind components to cell-center values with xarray. The destaggering is something that we can automatically handle with xwrf, but we explain the process below to be more transparent, especially with winds that are not in the cell-centers.

Destagger and Rename Dimensions

We need to:

  1. destagger to cell centers

  2. rename the staggered dimension back to the non-staggered name to avoid dimension conflicts

  3. (re)name the unstaggered wind for convenience

  4. Then, we are able to put these new DataArrays back into the ds_wrf Dataset.

# Plot wind vectors at a selected level to demonstrate how to destagger the wind components to cell-center values with xarray...

plot_level = 12  # index of level to plot
skip_xy = 10  # Sampling interval for the vector thinning

nt, nz, ny, nx = ds_wrf["T"].shape

ds_wrf["UA"] = 0.5*( ds_wrf["U"].isel(west_east_stag=slice(0, nx)) + 
                    ds_wrf["U"].shift(west_east_stag=-1).isel(west_east_stag=slice(0, nx)) ).\
                    rename("UA").rename(west_east_stag="west_east")
ds_wrf["VA"] = 0.5*( ds_wrf["V"].isel(south_north_stag=slice(0, ny)) + 
                    ds_wrf["V"].shift(south_north_stag=-1).isel(south_north_stag=slice(0, ny)) ).\
                    rename("VA").rename(south_north_stag="south_north")
ds_wrf["SPD"] = np.sqrt(ds_wrf["UA"]**2 + ds_wrf["VA"]**2).rename("wind speed").\
                    assign_attrs(units="m s-1", description="wind speed")

Visualize the Speed and Wind Vectors

Now, we can proceed to more plotting-specific data manipulation. We need to add spatial variables for the idealized domain (since XLAT and XLONG are constant in the file). This is needed by the xarray quiver routine.

Then, thin the grid to reduce the number of arrrows.

ds_wrf["west_east"] = xr.DataArray(data=np.arange(nx)*ds_wrf.attrs["DX"], dims="west_east", name="west_east", attrs={"units": "m"})
ds_wrf["south_north"] = xr.DataArray(data=np.arange(ny)*ds_wrf.attrs["DX"], dims="south_north", name="south_north", attrs={"units": "m"})

ds_wrf_thinned = ds_wrf.\
    isel(west_east=slice(0, nx, skip_xy), south_north=slice(0, ny, skip_xy), bottom_top=plot_level).\
    sel(Time=f"{case_date:%Y-%m-%d} {hour_to_plot}:00")

fig, ax = plt.subplots(ncols=1)
ds_wrf_thinned["SPD"].plot(ax=ax, x="west_east", y="south_north")
ds_wrf_thinned.plot.quiver(ax=ax, x="west_east", y="south_north", u="UA", v="VA",
                           scale=100)
plt.show()

Compare with Observational Data from ARM

Now that we have plotted the simulation data from WRF, let’s take a look at the observations.

fs.glob(f"{path}/*")
fs.glob('pythia/lasso-sgp/high_res_obs/sgp*')
# Compare with ARM Observational Data

import os

from arm_test_data import DATASETS
import matplotlib.pyplot as plt

import act

# Place your username and token here
username = os.getenv('ARM_USERNAME')
token = os.getenv('ARM_PASSWORD')

# If the username and token are not set, use the existing sample file
if username is None or token is None or len(username) == 0 or len(token) == 0:
    filename_ceil = DATASETS.fetch('sgpceilC1.b1.20190101.000000.nc')
    ceil_ds = act.io.arm.read_arm_netcdf(filename_ceil, engine='netcdf4')
else:
    # Example to show how easy it is to download ARM data if a username/token are set
    results = act.discovery.download_arm_data(
        username, token, 'sgpceilC1.b1', '2022-01-14', '2022-01-19'
    )
    ceil_ds = act.io.arm.read_arm_netcdf(results)

# Adjust ceilometer data for plotting
ceil_ds = act.corrections.ceil.correct_ceil(ceil_ds, -9999.0)

# Plot up ceilometer backscatter using HomeyerRainbow CVD friendly colormap
# The same could be done with keyword 'cmap='HomeyerRainbow'
display = act.plotting.TimeSeriesDisplay(ceil_ds, subplot_shape=(1,), figsize=(15, 5))
display.plot('backscatter', subplot_index=(0,), cvd_friendly=True)
plt.show()
# ARM Plotting v2
import act
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

Plot Desired Variables

# Set your username and token here!
username = '***'
token = '***'

# Set the datastream and start/enddates
datastream = 'sgpaosmetE13.a1'
startdate = '2019-04-04'
enddate = '2019-04-05'

# Use ACT to easily download the data.  Watch for the data citation!  Show some support
# for ARM's instrument experts and cite their data if you use it in a publication
result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)
# Let's read in the data using ACT and check out the data
ds_mpl = act.io.read_arm_netcdf(result)

ds_mpl

Change variable (as desired), Apply QC, and Plot Again

# Let's take a look at the quality control information associated with a variable from the MPL
variable = 'temperature_ambient'

# First, for many of the ACT QC features, we need to get the dataset more to CF standard and that
# involves cleaning up some of the attributes and ways that ARM has historically handled QC
ds_mpl.clean.cleanup()

# Apply corrections for the ceilometer, correcting for the vertical height
#ds_mpl = act.corrections.ceil.correct_ceil(ds_mpl,-999.0)

# Next, let's take a look at visualizing the quality control information
# Create a plotting display object with 2 plots
display = act.plotting.TimeSeriesDisplay(ds_mpl, figsize=(10, 5), subplot_shape=(1,))

# # Plot up the variable in the first plot
display.plot(variable, subplot_index=(0,))

# # Plot up the QC variable in the second plot
# #display.qc_flag_block_plot(variable, subplot_index=(1,))
# plt.show()
Project Pythia Cookbook Template
Project Pythia Cookbook Template
LASSO Simulation April 4, 2019 - Cloud base and Top
Clouds over SGP for April 4, 2019