Area averaged

Community matplotlib-users
1

Hi,
I just started using Python for the last few weeks, and previously been using GrADS for around 4 years.

I have trouble looking for a simplest way to calculate area average, let say I need to calculate a SST over a region of 20S-10N and 130E-170E.

I know how to get one point values of SST vs Time, as in:

···

…

…

f = nc.Dataset(‘d:/data/sst.mon.mean.nc’, ‘r’)

sstv = f.variables[‘sst’]

timev = f.variables[‘time’]

sst = sstv[:, 35, 100]

plt.plot(timev,sst)

plt.show()

…

…

but I couldn’t figure out how to get an area average value (…and didn’t get the right reference in the internet either)
Something missing, probably because I don’t understand enough about slicing or something else.

Can anyone give me a hint ?

Thanks.

Fadzil.

2

Hi Fadzil,

All you actually need to do is use numpy.average(), which is numpy’s implementation of the weighted average. It can be shown geometrically that using the cosine of the latitude as the weights in the weighted average would give you approximately the area average, though if your SST data has a grid cell area attribute in the netcdf file, that would be the most suitable choice to use as your weights. Otherwise, you could determine the area weighted average as follows:

numpy is imported as np, lat are the latitudes extracted from the netcdf file

First we need to convert the latitudes to radians

latr = np.deg2rad(lat)

Use the cosine of the converted latitudes as weights for the average

weights = np.cos(latr)

Assuming the shape of your data array is (nTimes, nLats, nLons)

First find the zonal mean SST by averaging along the latitude circles

sst = sstv[:]

sst_ave_zonal = sst.mean(axis=2)

Then take the weighted average of those using the weights we calculated earlier

sst_ave = np.average(sst_ave_zonal, axis=1, weights=weights)

This should give a time series of global mean SST. Is this what you wanted?

Thanks,

Alex

···

On Sun, Feb 23, 2014 at 10:28 AM, Fadzil Mnor <fadzilmnor84@…287…> wrote:

Hi,
I just started using Python for the last few weeks, and previously been using GrADS for around 4 years.

I have trouble looking for a simplest way to calculate area average, let say I need to calculate a SST over a region of 20S-10N and 130E-170E.

I know how to get one point values of SST vs Time, as in:

…

…

f = nc.Dataset(‘d:/data/sst.mon.mean.nc’, ‘r’)

sstv = f.variables[‘sst’]

timev = f.variables[‘time’]

sst = sstv[:, 35, 100]

plt.plot(timev,sst)

plt.show()

…

…

but I couldn’t figure out how to get an area average value (…and didn’t get the right reference in the internet either)
Something missing, probably because I don’t understand enough about slicing or something else.

Can anyone give me a hint ?

Thanks.

Fadzil.

Managing the Performance of Cloud-Based Applications

Take advantage of what the Cloud has to offer - Avoid Common Pitfalls.

Read the Whitepaper.

http://pubads.g.doubleclick.net/gampad/clk?id=121054471&iu=/4140/ostg.clktrk

Matplotlib-users mailing list

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

–
Alex Goodman
Graduate Research Assistant

Department of Atmospheric Science

Colorado State University

3

Thanks Alex for the reply.
So, that script calculates the global SST. What if when we want to calculate for only in specific box? For example, SST over this area only:

----------------------------------- 10 N

                                        >

                                        >

                                        >

              SST                  |

                                        >

                                        >

----------------------------------- 20 S

130 E 170E

Thanks.

···

Fadzil

4

Hi Fadzil,

Ah sorry, I glossed over that part of your question. There are actually two solutions to this, one would be to actually find the indices where the latitudes and longitudes are within your desired bounds using numpy.where(). However I generally prefer to use numpy’s built-in fancy indexing for this type of problem. For example:

lat and lon are extracted from the netcdf file, assumed to be 1D

Determine which latitudes are between 20S and 10N

latidx = (lat >= -20) & (lat <= 10)

Determine which longitudes are between 130E and 170E.

The numbers here may differ depending on the longitude convention in your data.

lonidx = (lon >= 130) & (lon <= 170)

Now we will actually subset the data. We need to subset lat too to make sure weights are consistent.

sst = sstv[:]

sst = sst[:, latidx][…, lonidx]

lat = lat[latidx]

···

Yes, the indexing does get a little tricky but it should work if you do it this way, then follow the same procedure outlined in the previous email.

Thanks,

Alex

On Sun, Feb 23, 2014 at 6:02 PM, Fadzil Mnor <fadzilmnor84@…287…> wrote:

Thanks Alex for the reply.
So, that script calculates the global SST. What if when we want to calculate for only in specific box? For example, SST over this area only:

----------------------------------- 10 N

| |

| |

| |

| SST |

| |

| |

----------------------------------- 20 S

130 E 170E

Thanks.

–
Alex Goodman
Graduate Research Assistant

Department of Atmospheric Science

Colorado State University

Fadzil