Reading multiple netCDF files and displaying monthly average data using Python

Monthly Average Sea Ice Concentration

Introduction

Here, we are going to use the Arctic sea ice concentration (SIC) data from the ASI (the Arctic Radiation And Turbulence Interaction Study (ARTIST) Sea Ice). The ASI algorithm computes daily sea ice concentration from the
polarization difference of 89 GHz brightness temperature. We are going to to calculate the monthly average of daily ASI SIC and display it using Python.

Data

The data files are in netCDF format (.nc). Daily ASI SIC product is published by the University of Bremen. The ASI SIC product has a polar stereographic projection at the 6.25 km spatial resolution and referenced to the WGS-84 ellipsoid, with each grid point center latitude and longitude provided with the data. We will be using data for the entire month of February, 2020. The ASI SIC netCDF files can be found here as “Feb 2020 ASI SIC netCDF Files”.

Necessary Libraries

NetCDF files can be read with a few different Python modules. We will use netCDF4. Follow this article for more about using netCDF library. We will also use numpy, maplotlib and glob libraries.

Reading netCDF file

We will use the netCDF.Dataset() attribute for reading the netCDF file. Let’s read a single netCDF file and explore the variables:

#Importing the necessary libraries
import numpy as np
import glob
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

#Reading a single netCDF file
data = Dataset(r'\path\to\file.nc')
print (data)

The output is:

Metadata of the netCDF file.

Here, ‘x( =1216)’, ‘y (=1792)’ and ‘z (= 0 to 100)’ variables are the number of columns, rows and the values of SIC, respectively.

First, let’s try to display a daily SIC data i.e., a single netCDF file.

#Creating a vairable 'sic' to store the 2-D array variable 'z'
sic = data.variables['z'][:]

#Displaying the variable 'z' i.e., SIC
cmmap = plt.cm.get_cmap("jet")
cmmap.set_bad('dimgrey',1.)
fig, ax = plt.subplots()
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.05)

im = ax.imshow(k, cmap=cmmap, vmin=0, vmax=100)

fig.colorbar(im, cax=cax, orientation='vertical', label='Sea Ice Concentration (%)')

ax.axes.xaxis.set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.gcf().set_size_inches(15, 15)
ax.set_title('ASI Daily SIC 1st Feb 2020 Grid 6.25 km', fontsize=20)
fig.set_dpi(300.0)
plt.show()

We get the output plot as:

We can see that the output SIC plot is flipped upside down. We can fix this using numpy.flip() to flip the image in an upright position. Let’s do this:

#Flipping the image upright in the axis = 0 i.e., vertically
f = np.flip(sic,0)

The resultant image is:

1st February ASI SIC.

Here, what we did was that we flipped the inverted image vertically i.e., in the axis = 0 position.

Now, let’s try to read the whole netCDF files for the entire month.

#File Path
w = glob.glob(r"path\to\files\*.nc")
#Creating an empty list to to append the SICs for each netCDF file
b = []
#Looping through the entire files and creating a variable for each SIC
for i in range(len(w)):
    exec (f'df{i} = Dataset(w[i])')
    exec (f"sea_ice_concentration{i} = df{i}.variables['z'][:]")
    exec (f"sic{i} = np.flip(sea_ice_concentration{i},0)")
    exec(f'b.append(sic{i})')

The *.nc in the file path is to read all the netCDF files. Notice the exec() function. It is useful for creating a new variable with a number at the end to keep track for each file. For example, the loop goes from 0 to 28 for this study and each netCDF file variable ‘z’ i.e., SIC values with 2-D arrays are created viz. sic0, sic1,sic2,sic3,…….. sic28. The same goes for ‘df’ and ‘sea_ice_concentration’ variables. The sic0, sic1,sic2,sic3,…….. sic28 variables so created are appended in the list b. The rationale behind this is to calculate the mean or average of all the daily SIC to get a monthly average SIC.

We are almost done now. We will use numpy.mean() to calculate the mean.

#Calculating monthly mean SIC
mean_sic = np.mean( np.array(b), axis=0 )

##Displaying the mean SIC
cmmap = plt.cm.get_cmap("jet")
cmmap.set_bad('dimgrey',1.)
fig, ax = plt.subplots()
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.05)

im = ax.imshow(mean_sic, cmap=cmmap, vmin=0, vmax=100)

fig.colorbar(im, cax=cax, orientation='vertical', label='Sea Ice Concentration (%)')

ax.axes.xaxis.set_visible(False)
ax.axes.yaxis.set_visible(False)
plt.gcf().set_size_inches(15, 15)
ax.set_title('ASI Monthly Average SIC Feb 2020 Grid 6.25km', fontsize=20)
#ax.add_artist(ScaleBar(10000))
fig.set_dpi(300.0)
plt.show()

In the above code, np.array(b) inside the np.mean() is to stacked the 2-D SICs above one another so that the corresponding pixels are aligned equally just like vertically stacking equally sized books upon one another.

The result:

Monthly Average ASI SIC February 2020

Conclusion

We have learnt how to plot the monthly average ASI SIC using python. This monthly average SIC data can be used for monthly average analysis of SIC. And eventually, we can use this for estimating sea ice thickness as SIC is one of the major factors for determining sea ice thickness.