Python-AWIPS Data Access

from awips.dataaccess import DataAccessLayer

# Connect to Unidata's cloud EDEX server
edex_url = "edex-cloud.unidata.ucar.edu"
DataAccessLayer.changeEDEXHost(edex_url)

print(f"Connected to EDEX server: {edex_url}")

# Get and display all available data types
data_types = DataAccessLayer.getSupportedDatatypes()
data_types.sort()
list(data_types)

# Create a new data request for grid data
request = DataAccessLayer.newDataRequest()
request.setDatatype("grid")

print(f"Created request for datatype: grid")

# Get available locations (models) for grid data
locations = DataAccessLayer.getAvailableLocationNames(request)
locations.sort()
list(locations)

# Set the location (model) we want to use
request.setLocationNames("GFS20")

print("Selected model: GFS20")

# Get available parameters for the selected model
params = DataAccessLayer.getAvailableParameters(request)
params.sort()

# Display first 30 parameters (there are many!)
print(f"Total parameters available: {len(params)}")
print("\nFirst 30 parameters:")
list(params[:30])

# Set the parameter we want
request.setParameters("T")

print("Selected parameter: T (Temperature)")

# Get available levels for the selected parameter
levels = DataAccessLayer.getAvailableLevels(request)

# Display first 20 levels
print(f"Total levels available: {len(levels)}")
print("\nSample levels:")
for lvl in levels[:20]:
    print(f"  {lvl}")

# Set the level we want
request.setLevels("0.0SFC")

print("Selected level: 0.0SFC (Surface)")

# Get available forecast cycles (run times)
cycles = DataAccessLayer.getAvailableTimes(request, True)

print(f"Number of available forecast cycles: {len(cycles)}")
print(f"\nMost recent cycle: {cycles[-1]}")

# Get all available times
times = DataAccessLayer.getAvailableTimes(request)

# Get the forecast run for the most recent cycle
fcst_run = DataAccessLayer.getForecastRun(cycles[-1], times)

print(f"Number of forecast times in latest run: {len(fcst_run)}")
print("\nFirst 10 forecast times:")
for t in fcst_run[:10]:
    print(f"  {t.getRefTime()} : {t.getFcstTime()} seconds ({t.getFcstTime()/3600:.0f} hours)")

# Request data for the first forecast time (analysis, 0 hour)
response = DataAccessLayer.getGridData(request, [fcst_run[0]])

print(f"Number of grid objects returned: {len(response)}")

# Extract information from the first grid object
grid = response[0]

print("Grid Metadata:")
print(f"  Model:     {grid.getLocationName()}")
print(f"  Parameter: {grid.getParameter()}")
print(f"  Level:     {grid.getLevel()}")
print(f"  Unit:      {grid.getUnit()}")
print(f"  Time:      {grid.getDataTime()}")
print(f"  Forecast:  {grid.getDataTime().getFcstTime()/3600:.0f} hours")

# Get the raw data array and coordinates
data = grid.getRawData()
lons, lats = grid.getLatLonCoords()

print(f"Data shape: {data.shape}")
print(f"Longitude shape: {lons.shape}")
print(f"Latitude shape: {lats.shape}")
print(f"\nData range: {data.min():.2f} to {data.max():.2f} {grid.getUnit()}")

import xarray as xr
import numpy as np

# Create an xarray DataArray from the grid data
da = xr.DataArray(
    data=data,
    dims=['y', 'x'],
    coords={
        'latitude': (['y', 'x'], lats),
        'longitude': (['y', 'x'], lons),
    },
    attrs={
        'units': grid.getUnit(),
        'long_name': grid.getParameter(),
        'level': str(grid.getLevel()),
        'model': grid.getLocationName(),
        'valid_time': str(grid.getDataTime()),
    },
    name=grid.getParameter()
)

da

from datetime import datetime, timedelta

# Request data for the first 5 forecast times
num_times = 5
response = DataAccessLayer.getGridData(request, fcst_run[:num_times])

# Build arrays for the Dataset
data_list = []
times_list = []

for grid in response:
    data_list.append(grid.getRawData())
    # Parse the time from the grid object
    ref_time = grid.getDataTime().getRefTime()
    fcst_seconds = grid.getDataTime().getFcstTime()
    # Create a proper datetime
    valid_time = datetime.strptime(str(ref_time)[:19], '%Y-%m-%d %H:%M:%S') + timedelta(seconds=fcst_seconds)
    times_list.append(valid_time)

# Stack data into a 3D array (time, y, x)
data_3d = np.stack(data_list, axis=0)

print(f"Combined data shape: {data_3d.shape}")

# Create the xarray Dataset
ds = xr.Dataset(
    data_vars={
        grid.getParameter(): (['time', 'y', 'x'], data_3d, {
            'units': grid.getUnit(),
            'long_name': 'Temperature',
            'level': str(grid.getLevel()),
        })
    },
    coords={
        'time': times_list,
        'latitude': (['y', 'x'], lats),
        'longitude': (['y', 'x'], lons),
    },
    attrs={
        'model': grid.getLocationName(),
        'source': 'EDEX via python-awips',
    }
)

ds

from metpy.units import units

ds = ds.metpy.parse_cf()

ds.T.metpy.units

# Select data by time
temp_first = ds['T'].isel(time=0)
print(f"Data at first time step:")
print(f"  Valid time: {temp_first.time.values}")
print(f"  Shape: {temp_first.shape}")

# Quick visualization
import matplotlib.pyplot as plt

# Convert from Kelvin to Celsius for plotting
temp_c = ds['T'].isel(time=0).metpy.convert_units('degC')

fig, ax = plt.subplots(figsize=(12, 8))
img = ax.pcolormesh(temp_c.longitude, temp_c.latitude, temp_c, cmap='RdBu_r', shading='auto')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title(f'GFS Surface Temperature\nValid: {temp_c.metpy.time.dt.strftime("%Y-%m-%d: %H:%MM:%SS")}')
cbar = plt.colorbar(img, ax=ax, label=f'Temperature [{temp_c.metpy.units:~P}]')
plt.show()