Xarray Fundamentals with Atmospheric Radiation Data

Xarray Fundamentals with Atmospheric Radiation Data#

This tutorial will introduce you to the Python library, Xarray. Xarray is a library designed for working with multi-dimensional arrays. This allows for the handling of geoscience data which typically involve multiple dimensions - for instance, latitude, longitude, time, and channels (Python Basics 5: Xarray — Digital Earth Africa Training 0.1 Documentation, n.d.). Xarray has two basic data types

xarray.DataArray - an array of data as well as its meta data

For example an array of temperatures across the surface of the earth.
xarray.Dataset - multiple data arrays merged into a single dataset with all associated metadata.

For example temperature and precipitation data across the surface of the earth.

Xarray is primary used to work with NetCDF files.

By the end of this tutorial, you would ahve gained a basic understanding of data manipulation and visualisation using Xarray.

We will use Xarray to analyze top-of-atmosphere radiation data from NASA’s CERES project.

Public domain, by NASA, from Wikimedia Commons

A pre-downloaded and subsetted a portion of the CERES dataset is available here: http://ldeo.columbia.edu/~rpa/CERES_EBAF-TOA_Edition4.0_200003-201701.condensed.nc. The size of the data file is 702.53 MB. It may take a few minutes to download.

Please review the CERES FAQs before getting started.

Start by importing Numpy, Matplotlib, and Xarray. Set the default figure size to (12, 6).

import xarray as xr
import numpy as np
from matplotlib import pyplot as plt
plt.rcParams['figure.figsize'] = (12, 6)
%config InlineBackend.figure_format = 'retina'

#filter warnings for cleaner notebook look
import warnings
warnings.filterwarnings('ignore')

Next, let’s download the NetCDF file.

Checking and changing working directory#

!pwd

Fetching the data#

import requests

data_url = "http://ldeo.columbia.edu/~rpa/CERES_EBAF-TOA_Edition4.0_200003-201701.condensed.nc"

response = requests.get(data_url, verify=False)

with open("CERES_EBAF-TOA_Edition4.0_200003-201701.condensed.nc", mode="wb") as file:
   file.write(response.content)

Code Explanation#

While pooch is a robust tool for retrieving data from remote sources, it does have limitations, particularly when facing issues such as errors in the URL hosting the dataset. In such cases, pooch might raise exceptions, necessitating an alternative approach.

To illustrate this, we define the URL of the data source and utilize the requests library to retrieve the data, setting SSL verification to False for cases where secure socket layer verification poses challenges.

data_url = "http://ldeo.columbia.edu/~rpa/CERES_EBAF-TOA_Edition4.0_200003-201701.condensed.nc"

response = requests.get(data_url, verify=False)

Subsequently, we save the file contents to a specified directory.

with open("CERES_EBAF-TOA_Edition4.0_200003-201701.condensed.nc", mode="wb") as file:
   file.write(response.content)

The data is now available with the specified name within the open() function.

To further illustrate this limitation, run the code below and compare the output to that using the alternative method.

import pooch
fname = pooch.retrieve(
    'http://ldeo.columbia.edu/~rpa/CERES_EBAF-TOA_Edition4.0_200003-201701.condensed.nc',
    known_hash='a876cc7106e7dcb1344fbec5dcd7510e5cd947e62049a8cbc188ad05ffe00345', 
    progressbar=True
)
print(fname)

Meridional Heat Transport Calculation#

Concepts for computation#

We can go beyond a weight factor and actually calculate the area of each pixel of the dataset. Because latitude and longitude correspond to spherical coordinates for Earth’s surface, each 1.0x1.0 degree grid cell actually has a different surface area as you move away from the equator! This is because latitudinal length is fixed (\( dLat = R d\varphi \)), but longitudinal length varies with latitude (\( d Lon = R d\lambda \cos(\varphi) \))

So the area element for lat-lon coordinates is

\[ dA = R^2 \cos(\lambda) d\lambda d \varphi \]

where \(\varphi\) is latitude, \(d\varphi\) is the spacing of the points in latitude, \(d\lambda\) is the spacing of the points in longitude, and \(R = 6,371\) km is Earth’s radius. (In this formula, \(\varphi\) and \(\lambda\) are measured in radians)

For more information, check out the xarray in 45 min by OceanHackWeek.

Calculating the pixel area using this formula and creating a 2D (lon, lat) DataArray for it#

We will need to specify the radius, \(d\varphi\) and \(d\lambda\) (from the dataset)

# Earth's average radius in km
R = 6371.0

# Coordinate spacing between each lat/lon is 1.0 degrees in the dataset
d_lambda = np.deg2rad(1.0)
d_phi = np.deg2rad(1.0)

#calculation
r_squared = R**2
dA = r_squared * np.cos(np.deg2rad(data.lat)) * d_lambda * d_phi

Now, to calculate the surface area

#broadcasting array
ones_array = xr.ones_like(data)
#multiply ones array by dA
pixel_area_da = ones_array * dA
#total_area
total_pixel_area = pixel_area_da.sum(dim=['lat', 'lon'])


#earth's surface area
earth_surface = 4 * np.pi * r_squared

The ones_like() function creates an array with the same shape and data type as the input array data, but with all elements set to 1.

Let’s compare our array calculation to the surface calculation using the surface area of a sphere formula.

#verify the calculation was correct

print(f" the pixel area is calculated to be {total_pixel_area.toa_sw_all_mon[0].values} and the surface area of the earth is {earth_surface}.")

 the pixel area is calculated to be 510070944.0 and the surface area of the earth is 510064471.90978825.

Calculating and plotting the total amount of net radiation in each 1-degree latitude band#

Multiplying the pixel area (m\(^2\)) from above with the radiative flux (W m\(^{-2}\)) gives you the total amount of radiation absorbed in each pixel in W.

#create dataset of area * net
net_by_area = dA * data.toa_net_all_mon


#create figure
plt.figure(figsize=(5,10))

#plot data - summing along lon and time leaving lat uncompressed
net_by_area.sum(dim=['lon','time']).plot(y='lat', color='deepskyblue',label='Net Radiation in Watts')

#major lines of latitude
plt.axhline(0, color='red', linestyle='solid', label='Equator')
plt.axhline(23.5, color='orangered', linestyle='dotted', label='Tropic of Cancer')
plt.axhline(-23.5, color='darkorange', linestyle='dotted', label='Tropic of Capricorn')
plt.axhline(66.5, color='cyan', linestyle='dashdot', label='Arctic Circle')
plt.axhline(-66.5, color='darkturquoise', linestyle='dashdot', label='Antarctic Circle')

#make legend and place outside of plot
plt.legend(loc='upper left', bbox_to_anchor=(1,1))

#create titles
plt.title('Total Absorbed Net Radiation')
plt.ylabel('Latitude')
plt.xlabel('Radiation in Watts (W)')

plt.grid(which='major', axis='x')

plt.show()

../../../_images/0b8523cbc7e25374d8f7a45ab5060e633ba6969a38c7cc1b6a1eb6b7fff051aa.png

Carefully review the code to generate this plot then make the same plot for incoming radation and total outgoing (short+long)

Plotting the cumulative sum of the total amount of net radiation as a function of latitude#

check out xarray’s cumsum function.

This curve tells you how much energy must be transported meridionally by the ocean and atmosphere in order to account for the radiative imbalance at the top of the atmosphere. Take a look at Chapter 2.1.5.2 of Introduction to Climate Dynamics and Climate Modelling for a better understanding of this concept.

A curve that looks something like this: https://journals.ametsoc.org/view/journals/clim/14/16/full-i1520-0442-14-16-3433-f07.gif (Figure from Trenberth & Caron, 2001)

# create a variable to hold the total
total = net_by_area.sum(dim=['lon','time'])

#take cumsum of total
cumsum_total = total.cumsum()

#plot cumsum
cumsum_total.plot()
#plot y=0
plt.axhline(0, color='red', linestyle='solid')

#set x limits
plt.xlim([-90,90])

#add grid
plt.grid(which='major', axis='y')

#create titles
plt.title('Cumulative Sum of Total Net Radiation by Latitude')
plt.xlabel('Degrees Latitude')
plt.ylabel('Heat Transport (PW)')

plt.show()

../../../_images/2e61b93eace06f31b689d7988d36df08011d007a24feb44f8f75cfffbba993a4.png

Selecting and Merging Data#

For the next problem, use the following approximate locations of four different cities.

city	lon	lat
Toronto, Ontario	79 W	42 N
Bridgetown, Barbados	59 W	13 N
Nome, Alaska	165 W	64 N
Columbo, Sri Lanka	80 E	7 N
Hobart, Tasmania	147 E	43 S

Remember the range of the index values for lon and lat, especially if slicing.

lon is \( [0, 360] \)
lat is \( [-90, 90]\)

Plot the incoming solar and net radiation at each city#

#index of the cities
index_list = list(Cities['City'].indexes['City'])

# Create subplots for each city

rows = len(index_list)
fig, axes = plt.subplots(rows, 1, figsize=(8, rows*4), sharex=False)


#iterate through the index list and select the city for plotting
for i, city in enumerate(index_list):
    data_net = Cities.toa_net_all_mon.sel(City=city)
    data_in = Cities.solar_mon.sel(City=city)

    # Plot each city's data on a separate subplot
    axes[i].plot(Cities.time, data_net, linewidth=0.8, label=f'{city} - TOA Net')
    axes[i].plot(Cities.time, data_in, linewidth=0.8, label=f'{city} - Incoming Solar')
    axes[i].axhline(0, linestyle='dashed', linewidth=0.3, color='black')
    axes[i].set_title(f"Timeseries of TOA Net and Solar Flux for {city}")
    axes[i].legend(bbox_to_anchor=(1,1.4))


plt.axhline(0, linestyle ='dashed', linewidth=0.3, color='black')
plt.xlabel('Time in Years')
#plt.legend(bbox_to_anchor=(1,1.4))
plt.tight_layout()
plt.show()

../../../_images/f7fd473afd5e41e6b11f4b1991cc7b2ea76b4c51b37a41335301259a2738436f.png

Examine the identified pattern and assess its suitability within the context of the specific location.

Final Thoughts#

By now, you should have acquired a foundational understanding of data structures and manipulation in Xarray through this example. If you have additional questions, I suggest referring to the comprehensive Xarray documentation for further guidance.

previous

Pandas Groupby with Hurricane Data

next

More Xarray with El Niño-Southern Oscillation (ENSO) Data

Contents

Checking and changing working directory

Fetching the data

Code Explanation

Opening data and examining metadata

Opening the dataset and displaying its contents

Printing metadata of each variable

Basic reductions, arithmetic, and plotting

Calculating the time-mean of the entire dataset

2D plot of the the time-mean Top of Atmosphere (TOA) Longwave, Shortwave, and Incoming Solar Radiation

Adding up the three variables above and verifying (visually) that they are equivalent to the TOA net flux

Deriving Net Flux Using Incoming and Outgoing Radiation

Plotting Net Flux Using TOA VAriable

Mean and weighted mean

Calculating the global (unweighted) mean of TOA net radiation

Creating a weight array proportional to \(\cos(\lambda)\)

Recalculating the global mean TOA net radiation with the weight factor

Now that you have a weight factor, we can verify that the TOA incoming solar, outgoing longwave, and outgoing shortwave approximately match up with infographic shown in the first cell of the tutorial

Meridional Heat Transport Calculation

Concepts for computation

Calculating the pixel area using this formula and creating a 2D (lon, lat) DataArray for it

Calculating and plotting the total amount of net radiation in each 1-degree latitude band

Plotting the cumulative sum of the total amount of net radiation as a function of latitude

Selecting and Merging Data

Creating a Dataset for each point from the global dataset

Merging these datasets into a new dataset with the new dimension city

Plot the incoming solar and net radiation at each city

Final Thoughts

By Marble Community

© Copyright 2024.

Xarray Fundamentals with Atmospheric Radiation Data

Contents

Xarray Fundamentals with Atmospheric Radiation Data#

Checking and changing working directory#

Fetching the data#

Code Explanation#

Opening data and examining metadata#

Opening the dataset and displaying its contents#

Printing metadata of each variable#

Basic reductions, arithmetic, and plotting#

Calculating the time-mean of the entire dataset#

2D plot of the the time-mean Top of Atmosphere (TOA) Longwave, Shortwave, and Incoming Solar Radiation#

Adding up the three variables above and verifying (visually) that they are equivalent to the TOA net flux#

Deriving Net Flux Using Incoming and Outgoing Radiation#

Plotting Net Flux Using TOA VAriable#

Mean and weighted mean#

Calculating the global (unweighted) mean of TOA net radiation#

Creating a `weight` array proportional to \(\cos(\lambda)\)#

Recalculating the global mean TOA net radiation with the weight factor#

Now that you have a `weight` factor, we can verify that the TOA incoming solar, outgoing longwave, and outgoing shortwave approximately match up with infographic shown in the first cell of the tutorial#

Meridional Heat Transport Calculation#

Concepts for computation#

Calculating the pixel area using this formula and creating a 2D (lon, lat) DataArray for it#

Calculating and plotting the total amount of net radiation in each 1-degree latitude band#

Plotting the cumulative sum of the total amount of net radiation as a function of latitude#

Selecting and Merging Data#

Creating a `Dataset` for each point from the global dataset#

Merging these datasets into a new dataset with the new dimension `city`#

Plot the incoming solar and net radiation at each city#

Final Thoughts#

Xarray Fundamentals with Atmospheric Radiation Data

Contents

Xarray Fundamentals with Atmospheric Radiation Data#

Checking and changing working directory#

Fetching the data#

Code Explanation#

Opening data and examining metadata#

Opening the dataset and displaying its contents#

Printing metadata of each variable#

Basic reductions, arithmetic, and plotting#

Calculating the time-mean of the entire dataset#

2D plot of the the time-mean Top of Atmosphere (TOA) Longwave, Shortwave, and Incoming Solar Radiation#

Adding up the three variables above and verifying (visually) that they are equivalent to the TOA net flux#

Deriving Net Flux Using Incoming and Outgoing Radiation#

Plotting Net Flux Using TOA VAriable#

Mean and weighted mean#

Calculating the global (unweighted) mean of TOA net radiation#

Creating a weight array proportional to \(\cos(\lambda)\)#

Recalculating the global mean TOA net radiation with the weight factor#

Now that you have a weight factor, we can verify that the TOA incoming solar, outgoing longwave, and outgoing shortwave approximately match up with infographic shown in the first cell of the tutorial#

Meridional Heat Transport Calculation#

Concepts for computation#

Calculating the pixel area using this formula and creating a 2D (lon, lat) DataArray for it#

Calculating and plotting the total amount of net radiation in each 1-degree latitude band#

Plotting the cumulative sum of the total amount of net radiation as a function of latitude#

Selecting and Merging Data#

Creating a Dataset for each point from the global dataset#

Merging these datasets into a new dataset with the new dimension city#

Plot the incoming solar and net radiation at each city#

Final Thoughts#

Creating a `weight` array proportional to \(\cos(\lambda)\)#

Now that you have a `weight` factor, we can verify that the TOA incoming solar, outgoing longwave, and outgoing shortwave approximately match up with infographic shown in the first cell of the tutorial#

Creating a `Dataset` for each point from the global dataset#

Merging these datasets into a new dataset with the new dimension `city`#