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.7487793This 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
Zoe's Insights 
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
LikeLike
I’m not a fan of Python. It seems to be popular among climate researchers and it’s faster than what I use, and so I learned a bit.
LikeLiked by 1 person
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?
LikeLike
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.
LikeLike