Skip to content

Commit

Permalink
ADD: New example showing how to use ACT to work with ARM and Satellit…
Browse files Browse the repository at this point in the history
…e data
  • Loading branch information
AdamTheisen committed Oct 11, 2024
1 parent 2a4ef57 commit f3266dc
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 14 deletions.
46 changes: 32 additions & 14 deletions act/plotting/xsectiondisplay.py
Original file line number Diff line number Diff line change
Expand Up @@ -146,13 +146,14 @@ def set_yrng(self, yrng, subplot_index=(0,)):

def plot_xsection(
self,
dsname,
varname,
field,
dsname=None,
x=None,
y=None,
subplot_index=(0,),
sel_kwargs=None,
isel_kwargs=None,
set_title=None,
**kwargs,
):
"""
Expand All @@ -162,11 +163,10 @@ def plot_xsection(
Parameters
----------
dsname : str or None
The name of the datastream to plot from. Set to None to have
ACT attempt to automatically detect this.
varname : str
field : str
The name of the variable to plot.
dsname : str or None
The name of the datastream to plot from.
x : str or None
The name of the x coordinate variable.
y : str or None
Expand All @@ -181,6 +181,8 @@ def plot_xsection(
to slice datasets, see the documentation on :func:`xarray.DataArray.sel`.
isel_kwargs : dict
The keyword arguments to pass into :py:func:`xarray.DataArray.sel`
set_title : str
Title for the plot
**kwargs : keyword arguments
Additional keyword arguments will be passed into
:func:`xarray.DataArray.plot`.
Expand Down Expand Up @@ -220,9 +222,9 @@ def plot_xsection(
y_coord_dim = temp_ds[y].dims[0]
coord_list[y] = y_coord_dim
new_ds = data_utils.assign_coordinates(temp_ds, coord_list)
my_dataarray = new_ds[varname]
my_dataarray = new_ds[field]
else:
my_dataarray = temp_ds[varname]
my_dataarray = temp_ds[field]

coord_keys = [key for key in my_dataarray.coords.keys()]
# X-array will sometimes shorten latitude and longitude variables
Expand Down Expand Up @@ -255,28 +257,42 @@ def plot_xsection(
self.set_xrng(xrng, subplot_index)
yrng = self.axes[subplot_index].get_ylim()
self.set_yrng(yrng, subplot_index)

if set_title is None:
if 'long_name' in self._ds[dsname][field].attrs:
set_title = self._ds[dsname][field].attrs['long_name']
plt.title(set_title)

del temp_ds
return self.axes[subplot_index]

def plot_xsection_map(
self, dsname, varname, subplot_index=(0,), coastlines=True, background=False, **kwargs
self,
field,
dsname=None,
subplot_index=(0,),
coastlines=True,
background=False,
set_title=None,
**kwargs,
):
"""
Plots a cross section of 2D data on a geographical map.
Parameters
----------
dsname : str or None
The name of the datastream to plot from. Set to None
to have ACT attempt to automatically detect this.
varname : str
field : str
The name of the variable to plot.
dsname : str or None
The name of the datastream to plot from.
subplot_index : tuple
The index of the subplot to plot inside.
coastlines : bool
Set to True to plot the coastlines.
background : bool
Set to True to plot a stock image background.
set_title : str
Title for the plot
**kwargs : keyword arguments
Additional keyword arguments will be passed into
:func:`act.plotting.XSectionDisplay.plot_xsection`
Expand All @@ -292,7 +308,9 @@ def plot_xsection_map(
'Cartopy needs to be installed in order to plot ' + 'cross sections on maps!'
)
self.set_subplot_to_map(subplot_index)
self.plot_xsection(dsname, varname, subplot_index=subplot_index, **kwargs)
self.plot_xsection(
field, dsname=dsname, subplot_index=subplot_index, set_title=set_title, **kwargs
)
xlims = self.xrng[subplot_index].flatten()
ylims = self.yrng[subplot_index].flatten()
self.axes[subplot_index].set_xticks(np.linspace(round(xlims[0], 0), round(xlims[1], 0), 10))
Expand Down
52 changes: 52 additions & 0 deletions examples/plotting/plot_satellite.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
"""
Using ACT for Satellite data
----------------------------
Simple example for working with satellite data.
It is recommended that users explore other
satellite specific libraries suchas SatPy.
Author: Adam Theisen
"""

import act
import matplotlib.pyplot as plt
import numpy as np
from arm_test_data import DATASETS

# Read in VISST Data
files = DATASETS.fetch('enavisstgridm11minnisX1.c1.20230307.000000.cdf')
ds = act.io.read_arm_netcdf(files)

# Plot up the VISST cloud percentage using XSectionDisplay
display = act.plotting.XSectionDisplay(ds, figsize=(10, 8))
display.plot_xsection_map(
'cloud_percentage', x='longitude', y='latitude', isel_kwargs={'time': 0, 'cld_type': 0}
)
plt.show()

# Download ARM TSI Data
files = DATASETS.fetch('enatsiskycoverC1.b1.20230307.082100.cdf')
ds_tsi = act.io.read_arm_netcdf(files)
ds_tsi = ds_tsi.where(ds_tsi.percent_opaque > 0)

# Set coordinates to extra data for ENA
ena_lat = 39.091600
ena_lon = 28.025700

lat = ds['lat'].values
lon = ds['lon'].values

# Find the nearest pixel for the satellite data and extract it
lat_ind = np.argmin(np.abs(lat - ena_lat))
lon_ind = np.argmin(np.abs(lon - ena_lon))

ds_new = ds.isel(lat=lat_ind, lon=lon_ind, cld_type=0)

# Plot the comparison using TimeSeriesDisplay
display = act.plotting.TimeSeriesDisplay({'Satellite': ds_new, 'ARM': ds_tsi}, figsize=(15, 8))
display.plot('cloud_percentage', dsname='Satellite', label='VISST Cloud Percentage')
display.plot('percent_opaque', dsname='ARM', label='ARM TSI Percent Opaque')
display.day_night_background(dsname='ARM')
plt.legend()
plt.show()

0 comments on commit f3266dc

Please sign in to comment.