20 June 2016

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

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

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



and surface:


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


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


  • ·       Downwelling longwave irradiance at surface

smaller ~16MB file:



larger (~50MB) one:



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:


Which is slightly different, but close enough, to the global TLT graph from RSS:



I may add more examples to this page as time goes on.


No comments:

Post a Comment