Planetary Vision: How to load .nc files into Matlab

Much climate data comes in the form of NetCDF files which have an .nc suffix.

This post is as much for my own reference in case I forget how to do it, as it is for anyone else interested in opening a NetCDF file in Matlab.

Here's a list of the current climate NetCDF files I am interested in:

· Satellite Temperature Lower Troposphere

http://www.remss.com/measurements/upper-air-temperature

↳

ftp://ftp.remss.com/msu/data/netcdf

↳

uat4_tb_v03r03_anom_chtlt_197812_201509.nc3.nc

(Note that the name suffix changes slightly with each month)

NCEP Reanalysis 2 of Upward longwave radiation flux at nominal top of atmosphere

http://iridl.ldeo.columbia.edu/SOURCES,/.NOAA/.NCEP-DOE/.Reanalysis-2/.Monthly/outline.html

http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCEP-DOE/.Reanalysis-2/.Monthly/.flx/.flx/.ulwrftoa/datafiles.html

and surface:

http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCEP-DOE/.Reanalysis-2/.Monthly/.flx/.flx/.ulwrfsfc/datafiles.html

· Satellite measurement of Outgoing Longwave radiation at nominal top of atmosphere.

http://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html

ftp://ftp.cdc.noaa.gov/Datasets/interp_OLR/olr.mon.mean.nc

· Downwelling longwave irradiance at surface

smaller ~16MB file:

http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.gaussian.html

http://www.esrl.noaa.gov/psd/repository/entry/show?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3AL25jZXAucmVhbmFseXNpczIuZGVyaXZlZC9nYXVzc2lhbl9ncmlkL2Rsd3JmLnNmYy5tb24ubWVhbi5uYw%3D%3D

larger (~50MB) one:

http://www.esrl.noaa.gov/psd/repository/entry/show?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3ALzIwdGhDX1JlYW5WMi9Nb250aGxpZXMvZ2F1c3NpYW4vbW9ub2xldmVsL2Rsd3JmLnNmYy5tb24ubWVhbi5uYw%3D%3D

With a few functions these files can be manipulated in Matlab such as: ncdisp, ncread, sum, NaNsum, and squeeze.

Let's start with RSS TLT. First we need to find the name of the variable to graph, which we can hopefully identify in the description that accompanies the data. We type

ncdisp('uat4_tb_v03r03_anom_chtlt_197812_201602.nc3.nc')

..which gives us the following output (I've deleted most lines for brevity):

Format:

classic

Global Attributes:

Conventions = 'CF-1.6'

title = 'Monthly anomalies of brightness temperatures in the Lower Troposphere on a 2.5 degree grid'

source = 'MSU and AMSU-A Level 1b'

references = 'doi:10.1175/2009JTECHA1237.1, doi:10.1175/2008JTECHA1176.1'

history = 'netcdf4 file converted from netcdf3 file written at Remote Sensing Systems'

.....

Variables:

longitude

Size: 144x1

Dimensions: longitude

Datatype: single

latitude

Size: 72x1

Dimensions: latitude

Datatype: single

time

Size: 458x1

Dimensions: time

Datatype: single

climatology_time

Size: 12x1

Dimensions: climatology_time

Datatype: single

longitude_bounds

Size: 2x144

Dimensions: bnds,longitude

Datatype: single

latitude_bounds

Size: 2x72

Dimensions: bnds,latitude

Datatype: single

time_bounds

Size: 2x458

Dimensions: bnds,time

Datatype: single

climatology_time_bounds

Size: 2x12

Dimensions: bnds,climatology_time

Datatype: single

brightness_temperature_anomaly

Size: 144x72x458

Dimensions: longitude,latitude,time

Datatype: single

brightness_temperature_climatology

Size: 144x72x12

Dimensions: longitude,latitude,climatology_time

Datatype: single

Notice brightness_temperature_anomaly has latitude, longitude and time in months, whereas brightness_temperature_climatology has only 12 months. The former is the quantity we want, while the latter is an average of all measured years over one 12 month period.

Now we can load the data in by defining a matrix, let's call it w1.

w1 =
ncread('uat4_tb_v03r03_anom_chtlt_197812_201602.nc3.nc','brightness_temperature_anomaly')

This yields a large (72x144x458) matrix in which we can view data for particular coordinates or a range of them.

Let's start with a single point at the South Pole. In the .nc file description under the variable latitude it says degrees north, where 1 is the South Pole and 72 is the North Pole.

So we want a small number for this southerly latitude. The satellite doesn't go all the way to latitude 1, so we have to start at about latitude = 9 to get some data.

For longitude it shouldn't matter which number we pick in the range 1 to 192, since on a Gaussian projection the points are close together at the poles, but let us choose long = 1, which is on the Greenwich meridian by definition.

w2 = w1(1,9,:)

This gives us the data we want but the matrix still has too many dimensions to graph so we use the squeeze function to reduce them:

w3 = squeeze(w2)

plot(w3)

This give us the following graph (which shows no warming in Antarctica since global warming™ began):

(X-axis is months from 1979, Y-axis is Temp anomaly (C))

Now let's try the North Pole and unlike Antarctica it should show some warming.

w4 = squeeze(w1(1,69,:))

Plot (w4)

And it shows the type of warming you'd expect at the North Pole. Warmists like to say that the North Pole is warming due to polar amplification, but Anarctica is betraying their script.)

This time instead of using a single point let's use all longitudes at the southern most latitude with data, latitude 9. First let's average out all longitudes to get rid of that dimension.

w5 = squeeze(sum(w1,1))

w5 is now a convenient matrix consisting of just bands of latitude.

Now we select the latitude 9:

w6 = w5(9,:)

..which gives us:

Now let's try a band of latitudes from 70S to 60S so we can compare it to the graph RSS gives.

w7 = sum(w5(9:13,:))

plot(w7)

This gives a good agreement with the RSS version (which has a few more months of data):

Here are the tropics from 20S to 20N. Note I have to use NaNsum this time because of some problem with the Not a Number entries.

w8 = nansum(w5(28:48,:))

And the RSS version:

The following is the band 60N to 82.5N near the North Pole. There's not quite as close an agreement to the official RSS graph:

w9 = nansum(w5(61:69,:))

So now that we have confirmed the coordinate system, at least in terms of latitude, let us try longitude.

The reason I mention this is that there is an error with one of these files – I can't remember which one - where the latitude is reversed. If I come across it I'll update the post, but I recommend confirming every coordinate system orientation if you can.

It wouldn't surprise me if errors like this end up in climate models. But anyway...

Let us dial in the continent of Australia. Using an Excel spread sheet I made as my Rosetta Stone to convert actual degrees of longitude and latitude (see here) to those in the NetCDF file, the range we need is: long 47 to 61 (113 East to 147 East) and lat 24 to 31 (17S to 34S)

w13 = squeeze(nansum(w1(47:61,:,:)))

w14 = nansum(w13(24:31,:))

plot(w14)

I haven't found suitable data to compare this to yet, but here's a similar graph from Steven Goddard:

A bit more here.

However the foregoing is not a perfect solution in that there are some minor differences in the global average I get when compared to the graphs officially produced.

To average across all lats and longs for a global average I use:

w11 = squeeze(nansum(w1,1))

w12 = squeeze(nansum(w11,1))

Here is the graph I get: