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.13/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.13/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.13/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.13/site-packages/s3fs/core.py:2378, in S3File._fetch_range(self, start, end)
   2377 try:
-> 2378     return _fetch_range(
   2379         self.fs,
   2380         self.bucket,
   2381         self.key,
   2382         self.version_id,
   2383         start,
   2384         end,
   2385         req_kw=self.req_kw,
   2386     )
   2388 except OSError as ex:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/site-packages/s3fs/core.py:2547, in _fetch_range(fs, bucket, key, version_id, start, end, req_kw)
   2546 logger.debug("Fetch: %s/%s, %s-%s", bucket, key, start, end)
-> 2547 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.13/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.13/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.13/site-packages/s3fs/core.py:2565, in _inner_fetch(fs, bucket, key, version_id, start, end, req_kw)
   2563         resp["Body"].close()
-> 2565 return await _error_wrapper(_call_and_read, retries=fs.retries)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/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.13/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.13/site-packages/s3fs/core.py:2552, in _inner_fetch.<locals>._call_and_read()
   2551 async def _call_and_read():
-> 2552     resp = await fs._call_s3(
   2553         "get_object",
   2554         Bucket=bucket,
   2555         Key=key,
   2556         Range="bytes=%i-%i" % (start, end - 1),
   2557         **version_id_kw(version_id),
   2558         **req_kw,
   2559     )
   2560     try:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/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.13/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.13/site-packages/xarray/backends/api.py:1788, in open_mfdataset(paths, chunks, concat_dim, compat, preprocess, engine, data_vars, coords, combine, parallel, join, attrs_file, combine_attrs, errors, **kwargs)
   1786 for p in paths1d:
   1787     try:
-> 1788         ds = open_(p, **open_kwargs)
   1789         datasets.append(ds)
   1790     except Exception as e:

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/site-packages/xarray/backends/api.py:760, 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)
    748 decoders = _resolve_decoders_kwargs(
    749     decode_cf,
    750     open_backend_dataset_parameters=backend.open_dataset_parameters,
   (...)    756     decode_coords=decode_coords,
    757 )
    759 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 760 backend_ds = backend.open_dataset(
    761     filename_or_obj,
    762     drop_variables=drop_variables,
    763     **decoders,
    764     **kwargs,
    765 )
    766 ds = _dataset_from_backend_dataset(
    767     backend_ds,
    768     filename_or_obj,
   (...)    779     **kwargs,
    780 )
    781 return ds

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/site-packages/xarray/backends/h5netcdf_.py:499, 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)
    496 emit_phony_dims_warning, phony_dims = _check_phony_dims(phony_dims)
    498 filename_or_obj = _normalize_filename_or_obj(filename_or_obj)
--> 499 store = H5NetCDFStore.open(
    500     filename_or_obj,
    501     format=format,
    502     group=group,
    503     lock=lock,
    504     invalid_netcdf=invalid_netcdf,
    505     phony_dims=phony_dims,
    506     decode_vlen_strings=decode_vlen_strings,
    507     driver=driver,
    508     driver_kwds=driver_kwds,
    509     storage_options=storage_options,
    510 )
    512 store_entrypoint = StoreBackendEntrypoint()
    514 ds = store_entrypoint.open_dataset(
    515     store,
    516     mask_and_scale=mask_and_scale,
   (...)    522     decode_timedelta=decode_timedelta,
    523 )

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

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/site-packages/xarray/core/utils.py:701, in read_magic_number_from_file(filename_or_obj, count)
    699 if filename_or_obj.tell() != 0:
    700     filename_or_obj.seek(0)
--> 701 magic_number = filename_or_obj.read(count)
    702 filename_or_obj.seek(0)
    703 return magic_number

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/site-packages/fsspec/spec.py:2111, in AbstractBufferedFile.read(self, length)
   2108 if length == 0:
   2109     # don't even bother calling fetch
   2110     return b""
-> 2111 out = self.cache._fetch(self.loc, self.loc + length)
   2113 logger.debug(
   2114     "%s read: %i - %i %s",
   2115     self,
   (...)   2118     self.cache._log_stats(),
   2119 )
   2120 self.loc += len(out)

File ~/micromamba/envs/lasso-those-clouds-cookbook-dev/lib/python3.13/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.13/site-packages/s3fs/core.py:2389, in S3File._fetch_range(self, start, end)
   2378     return _fetch_range(
   2379         self.fs,
   2380         self.bucket,
   (...)   2385         req_kw=self.req_kw,
   2386     )
   2388 except OSError as ex:
-> 2389     if ex.args[0] == errno.EINVAL and "pre-conditions" in ex.args[1]:
   2390         raise FileExpired(
   2391             filename=self.details["name"], e_tag=self.details.get("ETag")
   2392         ) from ex
   2393     else:

TypeError: argument of type 'NoneType' is not 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