Python NetCDF Latitude Area Weighted Global Average

The earth is not flat, and it’s not a sphere. If you’re doing data analysis on climate netcdf files with the Python programming language, you may need to figure out how to properly average gridded data over oblate spheroid Earth. Unfortunately there is few to no simple guides on how to do this properly. Searching google with “python netcdf latitude weighted area average” yielded pretty much nothing useful. Adding “tutorial” finally go me somewhere. It got me here: The Correct Way to Average the Globe .

The tutorial is okay, but the author does things in the most convoluted way possible. The amount of code made my eyes glaze over, and so I abandoned using Python for this task. But recently I had some time to think and I came up with a very short solution. It is in fact not very different from the Awk code I always used in many of my posts.

I hope this post will become the go-to guide for Python users on this question.

Here’s how it’s done:

# Zoe Phin, 2022/05/03
# File: berkland.py
# Run: python3 berkland.py
# Output: berk.png

import os
import requests as r
import xarray as x  
import numpy as n
from matplotlib import pyplot as m

if not os.path.exists('berk.nc'):
    open('berk.nc','wb').write(r.get('http://berkeleyearth.lbl.gov/auto/Global/Gridded/Complete_TAVG_LatLong1.nc').content)

d = x.open_dataset('berk.nc')['temperature']

a = 6378.137; e = 1-6356.752**2/a**2; r = n.pi/180

by_lat = (a*r)**2*(1-e)*n.cos(r*d.latitude)/(1-e*n.sin(r*d.latitude)**2)**2

w = d.weighted(by_lat).mean({'longitude','latitude'})
    
m.plot(w.time.values, w.values)
m.title('Land-Only Temperature Anomaly')
m.savefig('berk.png')

This is what we get:

Yes, this is correct.

We can easily do a sanity check. The sum area of all the grid cells should add up to the surface area of the Earth:

# File: area.py
import xarray as x  
import numpy as n

d = x.open_dataset('berk.nc')['temperature']

a = 6378.137; e = 1-6356.752**2/a**2; r = n.pi/180

by_lat = (a*r)**2*(1-e)*n.cos(r*d.latitude)/(1-e*n.sin(r*d.latitude)**2)**2

s=0
for l in by_lat.values:
    s += l*360
print(s)
$ python3 area.py

510072158.7487793

This is the correct area in km2 when using this technique.

Please note that this code works as is because the grid cell is 1×1 degrees. The way you modify it for different resolution is to change the by_lat line:

by_lat = (a*r)**2 ...      # 1 degree
by_lat = (a*r*0.5)**2 ...  # 0.5 degree
by_lat = (a*r*0.25)**2 ... # 0.25 degree
... 

That’s it.

Enjoy doing things the easy way.

-Zoe

Published by Zoe Phin

https://phzoe.com

4 thoughts on “Python NetCDF Latitude Area Weighted Global Average

  1. I don’t code any more, stopped at Honeywell DAP 516 Assembler and Fortran 4, written on KSR-33 BUT but my past boss and friend Bob Denny codes. He has used Python, whatever that is, in his telescope pointing and scheduling software business DC-3 Dreams. He posted this, which you might understand. I did not. He just got a Hubble Award so might be quite clever, but he never says so…

    https://www.facebook.com/photo/?fbid=10160341048619783&set=a.47021669782

    Like

  2. Back in Domebook days folks would cheat and divvy up an icosahedron with spherical trig, plot the resulting vertices in 3-space, multiply the vertical coordinate values by less than 1 and that yielded info for remapping the triangles onto a prolate spheroid and using the distance formula to sort out surface lengths. But nobody uses Dymaxion maps. The Gloege approach seems kinda like finding the area of the flat blueprint projection of a roof. Then instead of multiplying by the sec of pitch angle like a builder buying shingles, he accounts for the curve change using that weighing formula to tweak the way latitude spreads out. Question from non-coder: Where is the part that distinguishes between land-only and sea area… Land_and_Ocean_LatLong1.nc? And what about the huge difference between land/sea ratios north and south of the equator?

    Like

    1. I use a different file “Complete_TAVG_LatLong1.nc”. I’m averaging area-weighted anamolies not temperatures.

      The land and ocean file has a land_mask variable that marks 1 for land, 0 for ocean.

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: