A+ A A-

Nicolas Fauchereau

Last Updated on Tuesday, 31 July 2012 14:49

Hits: 108

Plotting: Maps

This tutorial:

File I/O: Reading ASCII data

Plotting: time-series

File I/O: Reading binary data

Plotting: Maps

 

In this tutorial, we will plot the TMI / TRMM data we read in Tutorial 3. For all plotting purposes I use the matplotlib package. The matplotlib package is very easy to learn for matlab users. Most of the methods in matplotlib use the same syntax as those in Matlab. For example, contour plots can be done with the contour method. Filled contour can be plotted with contourf and vector plots can be drawn with the quiver function, just like in Matlab. I strongly encourage you to visit the matplotlib page. In particular the Gallery page from matplotlib is very useful. In there you will find example of pretty much all the plots you are ever going to need.

To plot maps, we will use the Basemap package. Basemap is a matplotlib addon toolkit which makes it easy to project and map 2-dimensional data. Once you have installed basemap, you can use the same functions as in matplotlib (plot, quiver, contour, contourf) to plot your data on maps.

Ok. So let'd get started by loading in the TMI/TRMM data using our previously defined readTmi1D function. Then we will plot the SST.

from courseFunctions import readTmi1D
pname = ""
fname="tmi_20100909v4"
# We read the sst collected during the ascending path
lon,lat,sst = readTmi1D(pname+fname,0,1)

Let's import the Basemap package:

from mpl_toolkits.basemap import Basemap

We now need to decide which map projection we would like to use. There is a long list of projections available on the Basemap package which can be seen in the Basemap help.

help(Basemap)

The default value for the projection is "Cylindrical Equidistant", which actually does not change your lat and lon at all. The first thing to do is to define the map projection. For this we need to specify the longitude and latitude boundaries of our map.

# Region definition (Agulhas) minLon = 16
maxLon = 35
maxLat = -20
minLat = -43
# Base Map set-up with Mercator projection

m = Basemap(llcrnrlon=minLon, \
llcrnrlat=minLat, \
urcrnrlon=maxLon, \
urcrnrlat = maxLat, \
resolution = 'h', \
projection = 'merc', \
lon_0 = minLon+(maxLon-minLon)/2, \
lat_0 = minLat+(maxLat-minLat)/2)

My data will be plotted on a Mercator grid. I have defined my map boundary using the region definition. The resolution for the Map data is set to high. This gives better defined coastline, countries etc.... Let us draw a contour of the coastlines and fill in a color for the continent and lakes.

m.drawcoastlines()
m.fillcontinents(color='beige',lake_color='aqua')

Let us add some rivers in a thin blue line.

m.drawrivers(linewidth=0.2,color='b')

Let's define parallels and meridian lines every 5 degrees and plot them on the map.

# Import all functionality of numpy package from numpy import *
# Define parallels anf meridians
parallels = arange(floor(minLat),ceil(maxLat),5.0)
meridians = arange(floor(minLon),ceil(maxLon),5.0)

Now that we have set-up the map, we want to plot the SST data we have read from before. I use the pcolormesh function of matplotlib to do that. It is like pcolor except much faster (equivalent to pcolor with shading flat option in Matlab). To plot on the map I have just laid out, I must project my longitude and latitude arrays to my map projection.

mlon, mlat = m(*meshgrid(lon,lat))

pcolormesh does not handle the NaNs in my data so I must mask my sst using the ma (Masked Array) module of numpy.

sst_masked = ma.masked_where(isnan(sst),sst)
m.pcolormesh(mlon,mlat,sst_masked,vmin=10,vmax=30)

vmin and vmax specify the SST range I want my color axis to span. We can check that by displaying a colorbar().

colorbar()

By default, the plot colormap will be set to the rainbow colormap. There are a lot more colormaps available from matplotlib and basemap. Here are links to images of the matplotlib and Basemap color tables. We are now going to use a colormap from the Basemap package called GMT_wysiwyg.

from mpl_toolkits.basemap import cm
colormap = cm.GMT_wysiwyg
# Make masked values white
colormap.set_bad('w', 1.0)
m.pcolormesh(mlon,mlat,sst_masked,cmap = colormap,vmin=10,vmax=30)
title('SST plot') colorbar()

You can download the example code from: plotTmi.py file

Go to top of page

Last Updated on Monday, 06 August 2012 17:21

Hits: 309

File I/O: Reading binary data

This tutorial:

File I/O: Reading ASCII data

Plotting: time-series

File I/O: Reading binary data

Plotting: Maps

 

In this tutorial we will read a file of binary data collected from the tropical Rainfall Measuring Mission (TRMM) satellite. A description of the binary data can be found on the RSS / TMI site. The daily data is the most complicated so that what we are going to read. You should then find it easy to read 3-daily , weekly or monthly TMI/TRMM data based on this tutorial.

The file we are using here was downloaded from the SSMI FTP site.

From the site, we see that the daily binary data file consists of fourteen 0.25 x 0.25 degree grid (1440,320) byte maps. Seven ascending maps in the following order: Time (T), Sea Surface Temperature (S), 10-meter Surface Wind Speed using 11 GHz (Z), 10-meter Surface Wind Speed using 37 GHz (W), Atmospheric Water Vapor (V), Cloud Liquid Water (L), and Rain Rate (R), are followed by seven descending maps in the same order.

This means that the data will have the following dimensions (14,1440,320). If we want to separate the data into ascending and descending paths (which I want), then we can order the data so it has dimensions (2,7,1440,320), where (2) represents the number of satellite paths (ascensing or descending), (7) indicates the number of variables in the file (1440) is the length of the longitude and (320) is the length of the latitude. Before we even start reading the data, we can define our longitude and lontitude variables. For this we use the numpy arange method.

# Number of longitude and latitude points
nlon = 1440
nlat = 320
# Define lats and lons
lat = arange(-39.875,40,0.25)
lon = arange(-179.875,180,0.25)

Ok, so the first thing to do to read the data in is to define where our file is on our sytem and to open a file pointer just like when reading the ASCII data in the 1st tutorial. Of course, make sure you have unzipped your binary file first.

pname=""
fname="tmi_20100909v4"
fid = open(pname+fname, 'rb')

The method we use is called fromfile. There are other ways to do this but fromfile is supposed to be very fast so that is the one I have chosen. We load the binary data into a variable called data and reshape it straight away to have the data in the correct shape, rather than in 1 long data stream.

data = fromfile(file=fid, dtype=uint8).reshape((2,7,nlat,nlon))

The data has been read in so now we can close the file.

fid.close()

Now we want to convert the data into meaningful geophysical values using the scaling factors and offsets for each variable. We also want to discard all data that is bad such as land, contaminated data etc... On the RSS / TMI site, it says that all values above 250 can be flagged as bad. So we will mask all values above or equal to 250. The mkas variable has the same dimension as our data.

mask=data>250

To extract say the SST variable collected during the satellite descending path we will do the following:

sst = data[1,1,:,:]*scale_factor+offset

We then assign a NaN to all the bad data.

sst[where(mask[1,1,:,:])]=nan

You should now have read the SST for the descending path. You can test that everything looks OK by plotting it using the imshow method from pylab.

from pylab import imshow, colorbar, show
imshow(sst,origin='lower');
colorbar()

Although we've done the hardest bit, there still remains one thing to do. We need to re-order the longitude and the data to go from -180 to 180 rather than from 0 to 360 degrees East. To do this this I am going to create an empty array the same shape as my SST and then I am going to cut my SST variable in the middle, put the second half in front (from 0 to 180 degrees East) and put the first half (from 0 to -180) last. I use the numpy empty_like method to initiate my new SST variable. Here is the code:

# Re-order my SST so that it goes from longitude -180 to 180 rather than 0 to 360
# Flip dimension of variable along longitude
newSST = empty_like(SST)
newSST[:,0:nlon/2] = SST[:,nlon/2:nlon]
newSST[:,nlon/2:nlon] = SST[:,0:nlon/2]

Now my SST variable correctly spans the longitudes from -180 to 180 degrees East. That's it. We can clean it up and put it all together in one function called readTmi1D, which we add to our other functions in the courseFunctions.py file.

Go to top of page

Last Updated on Monday, 06 August 2012 17:23

Hits: 265

Plotting: time-series

This tutorial:

File I/O: Reading ASCII data

Plotting: time-series

File I/O: Reading binary data

Plotting: Maps

 

In this tutorial we are going to use the data we have just read from our ASCII file in File I/O: Reading ASCII data example and then plot it in a time-series.

For all plotting functionality I use pylab, which comes with the matplotlib package.

First let's load the data using the function we have created in the File I/O: Reading ASCII data example.

from courseFunctions import readCo2
# Path and file name definition
pname = ''
fname = 'CPT_CO2_dm_95_07.txt'
# Read in data
time, x, y = readCo2(pname+fname)

Using pylab we just plot data without any formatting.

from pylab import *
plot(time,x,'b')
plot(time,y,'r')
show()

This is all right but we would like to have our x axis formatted to indicate dates rather than just numbers. For this we use the plot_date method of matplotlib. We create a figure object and define the layout of our axis.

fig = figure(figsize =(12., 9.),facecolor='white', edgecolor='black')
ax = fig.add_axes([0.05, 0.5, 0.85, 0.35],axisbg='white')

We convert our time object to a matplotlib time object. I call the matplotlib time object plotDate.

from datetime import datetime, timedelta
st = matplotlib.dates.date2num(datetime(1950,1,1))
plotDate=matplotlib.dates.num2date(time+st)

I plot my unfiltered data as a thin grey line and my filtered data as a thick black line.

ax.plot_date(plotDate,x, 'gray',linewidth = 1.0)
ax.plot_date(plotDate,y, 'k',linewidth = 1.5)

I want my x-axis to start on the 1-Jan 1995 and to end on the 31st of December 2007. I define 2 matplotlib dates called start_time and end_time to specify the axis-limits. If you wanted to plot data for a shorter period, you could change the definition of start_time and end_time to something else.

start_time=matplotlib.dates.date2num(datetime(1995,1,1))
end_time=matplotlib.dates.date2num(datetime(2007,12,31))

I now set the x-axis limits to start_time and end_time.

ax.set_xlim( time_start, time_end)

Now the trick is to use the following methods from matplotlib.dates.

from matplotlib.dates import YearLocator, MonthLocator, DateFormatter

We want to have a reference for every year and every month on our plot.

years    = YearLocator()   # every year
months   = MonthLocator()  # every month

The following defines the format for how our time is displayed on the x-axis. In this case, it will be a 4-digit year.

dateFmt = DateFormatter('%Y')

We want labels on the x-axis for each year.

# format the ticks
ax.xaxis.set_minor_locator(months)
ax.xaxis.set_major_locator(years)
ax.xaxis.set_major_formatter(dateFmt)
ax.xaxis.grid(True,'minor',linewidth=1)
ax.xaxis.grid(True,'major',linewidth=1,linestyle='-')
xticks(size=9)
yticks(size=9)

We can now print out the plot using the FigureCanvasAgg method.

from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
canvas = FigureCanvas(fig)
canvas.print_figure('co2_series.png',dpi=100)

You can download the code from file plotCo2.py.

Go to top of page

Last Updated on Monday, 06 August 2012 17:19

Hits: 215

CHPC Student Cluster Competition 2013

Tsessebe Cluster Available

Graphical Processing Unit Cluster Available

CHPC SAGrid Cluster Available

Dirisa Storage Unit Available

Social shares

Website developed by Multidimensions