Processing: CARRA Reanalysis Atmosphere

In this notebook, we’ll work with CARRA (Copernicus Arctic Regional Reanalysis) data to construct atmospheric profiles for pyRadtran simulations. CARRA provides higher-resolution atmospheric data for the Arctic compared to ERA5.

Note

Campaign data required. This notebook uses CARRA reanalysis data files that may not be publicly available in the specific format shown. You can adapt the workflow to your own CARRA or ERA5 data.

import pyradtran
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import logging
from pathlib import Path
from pyradtran.config import load_config

logging.getLogger('pyradtran').setLevel(logging.DEBUG)

Setup

kiruna_lat, kiruna_lon = 67.84889, 20.3027

tgt_store = "/projekt_agmwend/data_raw/COMPEX-EC/04_ERA5/CARRA_pressure_levels/CARRA_point_optimized.zarr"

Loading CARRA Data

ds_kiruna_regridded = xr.open_dataset('/projekt_agmwend/data_raw/COMPEX-EC/04_ERA5/CARRA_pressure_levels/CARRA_EAST_latlon0p05.nc', chunks='auto').rename({'lat' : 'latitude', 'lon': 'longitude'})

ds_kiruna_optimized = xr.open_dataset(tgt_store)
ds_kiruna = ds_kiruna_optimized.sel(lat=kiruna_lat, lon=kiruna_lon, method='nearest').load()
ds_kiruna.r.units

ds_kiruna_to_sim = ds_kiruna
ds_kiruna_to_sim = ds_kiruna_to_sim.rename({'r': 'q', 
                                            'plev': 'pressure_level',
                                            'lon' : 'longitude',
                                            'lat' : 'latitude',
                                            'time' : 'valid_time'})  # pyradtran expects 'q' for specific humidity
ds_kiruna_to_sim['pressure_level'].attrs['units'] = 'Pa'  # change unit from Pa to hPa

ds_kiruna_to_sim.pressure_level.units

from pyradtran.io import ERA5AtmosphereGenerator

atmosphere = ERA5AtmosphereGenerator()
atmosphere.create_era5_atmosphere_file(
    ds_kiruna_to_sim,
    latitude=kiruna_lat,
    
    longitude=kiruna_lon,
    time=ds_kiruna_to_sim.valid_time.values[0],
    output_filepath='work/test.atm'
)

ds_kiruna.clwc.plot(x=‘time’, yincrease=False) plt.yscale(‘log’)

import pandas as pd

# Build thermal configuration from master config
cfg = load_config()
cfg.simulation_defaults.rte_solver = "disort"
cfg.simulation_defaults.mol_abs_param = "reptran medium"
cfg.simulation_defaults.source = "thermal"
cfg.simulation_defaults.wavelength_nm = [3000, 100000]
cfg.simulation_defaults.integrate_wavelength = True
cfg.simulation_defaults.albedo_value = 0.85
cfg.simulation_defaults.surface_temperature_k = 248.4
cfg.simulation_defaults.output_altitudes_km = list(range(1, 11))
cfg.simulation_defaults.output_columns = ["zout", "lambda", "sza", "edir", "eglo", "edn", "eup", "enet", "albedo"]
cfg.execution.max_workers = 16
cfg.execution.cleanup_temp_files = False
cfg.output.filename_prefix = "pyradtran_sim"
config_path = Path("config/thermal_config.yaml")
cfg.to_yaml(config_path)

logging.getLogger('pyradtran').setLevel(logging.DEBUG)

# Stationary simulation at constant altitude, latitude, and longitude but with varying time
N_timesteps = ds_kiruna.dims['time']

ds = xr.Dataset(
    coords={
        'time' : ds_kiruna['time'],
        'latitude' : ('time', [kiruna_lat] * N_timesteps),
        'longitude' : ('time', [kiruna_lon] * N_timesteps),
        'altitude' : ('altitude', np.arange(0, 20, 1)),
    }
)

Running Simulations

ds_kiruna_to_sim

ds_sim = ds.pyradtran.run(
    config_path=config_path,
    return_dataset=True,
    save_to_file=False,
    era5_atmosphere=ds_kiruna_to_sim,  # Use the local ERA5 data
    show_progress=False,
)

Distance between Requested Point and Nearest ERA5 Point

Compare the distance between the requested point and the nearest ERA5 point with the haversine function.

def haversine(lon1, lat1, lon2, lat2):
    # convert decimal degrees to radians
    lon1 = np.deg2rad(lon1)
    lon2 = np.deg2rad(lon2)
    lat1 = np.deg2rad(lat1)
    lat2 = np.deg2rad(lat2)

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    r = 6371
    return c * r

distances = haversine(
    ds['longitude'].values,
    ds['latitude'].values,
    ds_kiruna_to_sim['longitude'].values,
    ds_kiruna_to_sim['latitude'].values
)

print(f"Distances to nearest ERA5 point: {distances.mean():.2f} km")

Results

As the simulation results are stored in an xarray dataset, you can use the xarray plotting capabilities to visualize the results.

ds_kiruna_to_sim

ds_kiruna_to_sim['z']

ds_kiruna_to_sim.pressure_level.units

ds_kiruna_to_sim['altitude']

fig, ax = plt.subplots(2, 1, figsize=(20, 6))

#change the prop cycle
ax[0].set_prop_cycle(plt.cycler(color=plt.cm.viridis(np.linspace(0, 1, N_timesteps))))
ax[1].set_prop_cycle(plt.cycler(color=plt.cm.viridis(np.linspace(0, 1, N_timesteps))))

ds_sim.enet.plot(
    x='time',
    y='altitude',
    ax=ax[1],
    cmap='viridis',
    robust=True

)

ds_kiruna_to_sim['altitude'] = ds_kiruna['z'] / 9.81 / 1e3
ds_kiruna_to_sim.assign_coords(altitude=ds_kiruna_to_sim['z'])

ds_kiruna_to_sim.t.plot.pcolormesh(
    x='valid_time',
    ax=ax[0],
    y='pressure_level',
    cmap='inferno',
    robust=True,
    yincrease=False
)
ax[0].set_yscale('log')