autoload.py

# coding: utf-8
import pprint
import yt
import numpy as np
import re

import netCDF4 as nc4

poor_regex = re.compile(r'([a-zA-Z]+)(.*)')
def poor_parser(s):
    try:
        return "**".join(filter(None, poor_regex.search(s).groups()))
    except AttributeError:
        return s

n = nc4.Dataset("./GEOS.fp.asm.inst3_3d_aer_Nv.20151030_0000.V01.nc4")

dims = []
sizes = []
bbox = []
ndims = len(n.dimensions)
for dim in n.dimensions.keys():
    size = n.variables[dim].size
    if size > 1:
        bbox.append([n.variables[dim][:].min(),
                     n.variables[dim][:].max()])
        dims.append(n.variables[dim].long_name)
        sizes.append(size)

dims.reverse()   # Fortran ordering
sizes.reverse()
bbox.reverse()
dims = [f.replace('vertical level', 'altitude') for f in dims]
bbox = np.array(bbox)

data = {}
names = {}
for field, d in n.variables.items():
    if d.ndim != ndims:
        continue
    units = n.variables[field].units
    units =  " * ".join(map(poor_parser, units.split())) 
    data[field] = (np.squeeze(d), str(units))
    names[field] = n.variables[field].long_name.replace("_", " ")

ds = yt.load_uniform_grid(data, sizes, 1.0, geometry=("geographic", dims),
                          bbox=bbox)
ds.index
for name in names:
    print ds.field_info["gas", name].display_name, "=>", names[name]
    ds.field_info["gas", name].display_name = names[name]
#ds.surface_height = 0 * ds.domain_width.uq
pprint.pprint(ds.field_list)
p = yt.SlicePlot(ds, "altitude", ds.field_list)
p.save("frames/")