#!/usr/bin/env python3 """The five on the measured planet, by ten-degree box. A place is counted twice. Once as ground: what the surface is made of, a hundred parts divided between earth and water, reckoned for a column going down from the surface to a fixed depth. Land is earth all the way down, less the water in it (see wetness below). The sea is water down to its floor, and below the floor it is soaked ground, half water and half earth; a sea deeper than the column is all water. Ice is water lying on soaked ground. So the shallowest sea and the wettest land meet at half and half. Once as light: where the sunlight that reaches it goes, of which three parts are kept, the part thrown back (space), the part handed to the air as warmth (air) and the part radiated away as heat (fire). Air is counted by its movement as well: half by that warmth and half by the wind, the wind set against the planet's mean wind and given the planet's mean warming for its unit. The parts of the sunlight that go into the ground itself, as evaporation or as stored heat, are not counted a second time. The light is counted by how much of it there is: each of the three, in watts, against the sunlight the planet's surface receives on average. So the ground of every place counts for a hundred, and the light of a place for more or less than the planet's forty-six according to how much sun it has. The two counts are added and the total is called a hundred. Reads ERA-Interim monthly surface fluxes and GPCP rainfall (and CERES-EBAF surface radiation as a check) from the CMIP Rapid Evaluation Framework's sample data, which carries them already regridded to points ten degrees apart, and the ground raster made by build_ground.py, and writes obs.json the shares themselves (for the page and for the game's starting map) obs_measured.json every number the two essays quote from the measurements Usage: python3 build_obs.py [--fetch] DATA_DIR OUT_DIR [ground_now.json] --fetch downloads the source files into DATA_DIR first (about 15 MB). Needs numpy and xarray (with netCDF4 or h5netcdf). Some things are tried beside the count and used nowhere in it (see "two tests" in main): a fire counted by how warm the surface is, a soil with a store that carries water from month to month, and the air counted by either of its halves alone. Two rules run through everything below. Over a region the light is pooled by sunlight (fluxes summed, then divided) and the ground by area, and the two are then joined as they are at a single box. Doshas are weaves, and a weave of a pool is not the pool of the weaves. Over a region the weave is made at each box from what is there, and the boxes are averaged by area. The one exception is the planet taken whole, which is given both ways. """ import base64, json, os, sys, urllib.request, warnings, zlib import numpy as np import xarray as xr warnings.filterwarnings('ignore') REPO = 'https://raw.githubusercontent.com/Climate-REF/ref-sample-data/08a1e3a209ad10acbb82a50c7b71822a1999a741/data/' ERA = 'obs4REF/obs4REF/ECMWF/ERA-INT/mon/{v}/gn/v20210727/{v}_mon_ERA-INT_PCMDI_gn_197901-201903.nc' CERES = 'obs4REF/obs4REF/NASA-LaRC/CERES-EBAF-4-2/mon/{v}/gn/v20230209/{v}_mon_CERES-EBAF-4-2_RSS_gn_200003-202309.nc' ERA_VARS = ['rsds', 'rsus', 'rlds', 'rlus', 'hfls', 'hfss'] ERA_EXTRA = ['ts', 'sfcWind'] # surface temperature (for comparison and one trial) and wind speed (half of the air) GPCP = 'obs4REF/obs4REF/NOAA-NCEI/GPCP-2-3/mon/pr/gn/v20210727/pr_mon_GPCP-2-3_PCMDI_gn_197901-201907.nc' CERES_VARS = ['rsds', 'rsus', 'rlds', 'rlus'] NAMES = ['space', 'air', 'fire', 'water', 'earth'] Y0, Y1 = 1979, 2018 # forty whole years NOW, THEN = (2009, 2018), (1979, 1988) # "today" and the first decade of the record DARK = 20 # W m-2: a month with less sunlight than this is counted as polar night (one statistic only) LV = 2.45e6 # J kg-1: the heat it takes to turn water to vapour GROUND = 1.0 # how many hundreds the ground counts for, against the light's one SOAK = 0.5 # soaked ground is half water: land is never less than half earth MOVING = 0.5 # how much of the air is counted by its movement (the wind); the rest by the heat handed to it COLUMN = 3700.0 # m: how far down the ground is counted; the mean depth of the ocean ICE = 2000.0 # m: the thickness taken for the ice sheets, which is near the mean of both # Share of the first strand: vata = air with space, pitta = fire with water, kapha = water with earth. WEAVE = dict(vata=3 / 4, pitta=1 / 2, kapha=2 / 3) TASTES = dict(vata=3 / 4, pitta=2 / 3, kapha=5 / 7) # the same count with the tastes of the lists included EQUAL = dict(vata=1 / 2, pitta=1 / 2, kapha=1 / 2) def fetch(data_dir): os.makedirs(data_dir, exist_ok=True) for tmpl, vs in ((ERA, ERA_VARS + ERA_EXTRA), (CERES, CERES_VARS), (GPCP, ['pr'])): for v in vs: p = tmpl.format(v=v) dst = os.path.join(data_dir, os.path.basename(p)) if not os.path.exists(dst): print('fetching', os.path.basename(p)) urllib.request.urlretrieve(REPO + p, dst) def load(data_dir, tmpl, v, a, b): ds = xr.open_dataset(os.path.join(data_dir, os.path.basename(tmpl.format(v=v)))) return ds[v].sel(time=slice(a, b)).values.astype('f8'), ds def pos(a): return np.clip(a, 0, None) def split(F): """Fluxes -> where the sunlight at the ground goes first (W m-2 each, summing to rsds). Works on whatever step the fluxes are given at: months here, whole years for the comparison.""" rsds = pos(F['rsds']) rsus = pos(np.minimum(F['rsus'], rsds)) A = rsds - rsus # absorbed by ground or sea LW = F['rlus'] - F['rlds'] # net thermal radiation, upward LH, SH = F['hfls'], F['hfss'] # evaporation; heat handed to the air G = A - LW - LH - SH # what is left goes into the ground or the sea out = pos(LW) + pos(LH) + pos(SH) + pos(G) f = np.clip(np.where(out > 0, A / np.where(out > 0, out, 1), 0), 0, 1) # the part the sun itself is paying for return dict(space=rsus, air=pos(SH) * f, fire=pos(LW) * f, evap=pos(LH) * f, stored=pos(G) * f, rsds=rsds, dry=pos(A - LW) / LV, # the sun's power to dry: what the ground absorbs and does not radiate, as water it could evaporate (kg m-2 s-1) raw=dict(rsds=rsds, rsus=rsus, lw=LW, lh=LH, sh=SH, g=G)) def light3(Pp): return np.stack([Pp['space'], Pp['air'], Pp['fire']], -1) def wetness(rain, dry, soak=None): """How much of a piece of land is water. Ground that gets as much rain as the sun could dry is soaked, and soaked ground is half water; ground that gets less is wet in proportion.""" soak = SOAK if soak is None else soak full = np.where(dry > 0, np.minimum(1.0, rain / np.where(dry > 0, dry, 1)), np.where(rain > 0, 1.0, 0.0)) return soak * full def join(l3, sun, water, earth, G=None): """Light (three fluxes and the sunlight they are shares of) and ground (two shares summing to one) -> the five.""" G = GROUND if G is None else G sun = np.asarray(sun, float) l = np.asarray(l3) / np.where(sun > 0, sun, 1)[..., None] l = np.where((sun > 0)[..., None], l, 0) if G == 'fixed': # the alternative rule: two-fifths to the ground, three-fifths to the light lt = l.sum(-1, keepdims=True) return np.concatenate([0.6 * l / np.where(lt > 0, lt, 1), 0.4 * np.stack([water, earth], -1)], -1) tot = l.sum(-1) + G return np.concatenate([l, G * np.stack([water, earth], -1)], -1) / tot[..., None] def weave(a, b, w): a, b = np.maximum(a, 0), np.maximum(b, 0) return (a / w) ** w * (b / (1 - w)) ** (1 - w) def strengths(E, W=WEAVE): """E[..., 5] (space, air, fire, water, earth) -> [..., 3] (vata, pitta, kapha), not yet set against one another.""" E = np.asarray(E) return np.stack([weave(E[..., 1], E[..., 0], W['vata']), weave(E[..., 2], E[..., 3], W['pitta']), weave(E[..., 3], E[..., 4], W['kapha'])], -1) def norm(S): S = np.asarray(S) t = S.sum(-1, keepdims=True) return np.where(t > 0, S / np.where(t > 0, t, 1), 1 / 3) def r(x, n=1): if isinstance(x, (list, tuple, np.ndarray)): return [r(v, n) for v in x] return round(float(x), n) def pct(x, n=1): if isinstance(x, (list, tuple, np.ndarray)): return [pct(v, n) for v in x] return round(float(x) * 100, n) def b64u16(a, scale=10000): return base64.b64encode(np.clip(np.round(a * scale), 0, 65535).astype(' 2 else os.path.join(os.path.dirname(os.path.abspath(__file__)), 'ground_now.json') F = {} for v in ERA_VARS: F[v], ds = load(data_dir, ERA, v, f'{Y0}-01', f'{Y1}-12') lat, lon = ds['lat'].values, ds['lon'].values NY, NLAT, NLON = Y1 - Y0 + 1, len(lat), len(lon) assert F['rsds'].shape == (NY * 12, NLAT, NLON) w = np.sin(np.radians(np.clip(lat + 5, -90, 90))) - np.sin(np.radians(np.clip(lat - 5, -90, 90))) W = np.repeat(w[:, None], NLON, 1) W /= W.sum() def gm(a): return float((a * W).sum()) P = split(F) L3 = light3(P) # month, lat, lon, 3 (W m-2, sun-paid) # Air is counted twice: by the heat the ground hands to it, and by how fast it moves. The wind is set against the planet's # mean wind and given the planet's mean heating for its unit, so that over the planet the two counts come to the same total. windM = load(data_dir, ERA, 'sfcWind', f'{Y0}-01', f'{Y1}-12')[0] # month, lat, lon, m s-1, ten metres above the ground n0_, n1_ = (NOW[0] - Y0) * 12, (NOW[1] - Y0 + 1) * 12 WBAR, ABAR = gm(windM[n0_:n1_].mean(0)), gm(L3[n0_:n1_, ..., 1].mean(0)) # the planet's mean wind, and its mean heat handed to the air, today HEAT, MOVE = L3[..., 1].copy(), ABAR * windM / WBAR def air_of(heat, move, m=None): m = MOVING if m is None else m return (1 - m) * heat + m * move L3[..., 1] = air_of(HEAT, MOVE) sun = P['rsds'] L3Y = L3.reshape(NY, 12, NLAT, NLON, 3) sunY = sun.reshape(NY, 12, NLAT, NLON) rainM = pos(load(data_dir, GPCP, 'pr', f'{Y0}-01', f'{Y1}-12')[0]) # kg m-2 s-1 rainY = rainM.reshape(NY, 12, NLAT, NLON) dryY = P['dry'].reshape(NY, 12, NLAT, NLON) # ---------- what the ground is ---------- d = json.load(open(ground_file)) res = d['res'] unpack = lambda k: np.unpackbits(np.frombuffer(base64.b64decode(d[k]), dtype=np.uint8))[:d['ny'] * d['nx']].reshape(d['ny'], d['nx']) rasters = {k: unpack(k) for k in ('land', 'ice', 'lake')} depth_class = np.frombuffer(zlib.decompress(base64.b64decode(d['depth'])), dtype=np.uint8).reshape(d['ny'], d['nx']) edges = [0] + list(d['depthLevels']) + [d['depthLevels'][-1] + 1000] # the depth classes, in metres for k in range(len(edges) - 1): rasters['sea%d' % k] = ((depth_class == k) & (rasters['land'] == 0)).astype(np.uint8) def at(k, la, lo): lo = ((lo + 180) % 360) - 180 return rasters[k][min(d['ny'] - 1, max(0, int((90 - la) / res))), int((lo + 180) / res) % d['nx']] def near_share(k, la, lo, half, n): # the share of the ground within `half` degrees of the point that is of kind k offs = np.linspace(-half, half, n) return np.mean([at(k, min(89.9, max(-89.9, la + dy)), lo + dx) for dy in offs for dx in offs]) # a point's fluxes were read from a degree and a half around it, so its ground is read from the same patch near = np.array([[near_share('land', la, lo, 1.5, 5) for lo in lon] for la in lat]) ice = np.array([[near_share('ice', la, lo, 1.5, 5) for lo in lon] for la in lat]) lake = np.array([[near_share('lake', la, lo, 1.5, 5) for lo in lon] for la in lat]) landfrac = np.array([[near_share('land', la, lo, 4.75, 20) for lo in lon] for la in lat]) # the whole ten-degree box poles = np.abs(lat) == 90 # a pole is one place, whatever longitude it is given for arr in (near, ice, lake, landfrac): arr[poles] = arr[poles].mean(1, keepdims=True) soil = np.clip(near - ice - lake, 0, 1) # bare or planted ground: earth, less the water in it sea_by_depth = np.array([[[near_share('sea%d' % k, la, lo, 1.5, 5) for lo in lon] for la in lat] for k in range(len(edges) - 1)]) sea_by_depth[:, poles] = sea_by_depth[:, poles].mean(2, keepdims=True) sea_by_depth[0] += 1 - soil - ice - lake - sea_by_depth.sum(0) # whatever rounding leaves over is shallow sea sea_by_depth[0] += lake # lakes are counted with the shallowest water def floor_share(column): """For each depth class, how much of a column of the given depth is sea floor: the mean of max(0, 1 - depth/column) across the class.""" out = [] if column <= 0: # the surface alone: no floor anywhere return np.zeros(len(edges) - 1) for lo_, hi_ in zip(edges[:-1], edges[1:]): dd = np.linspace(lo_, hi_, 201) out.append(float(np.clip(1 - dd / column, 0, None).mean())) return np.array(out) def under(column): """The ground under the water of each box, as a share of its column: sea floor within the column, and the land under the ice.""" return (sea_by_depth * floor_share(column)[:, None, None]).sum(0) + ice * (max(0.0, 1 - ICE / column) if column > 0 else 0.0) def ground(rain, dry, soak=None, column=None, wet=None): """Water and earth of each box as shares of its ground, for a period's summed rain and drying power. The ground under water is soaked, so only the part of it that is not water counts as earth. `wet`, if given, is the wetness of the soil reckoned some other way (the soil that remembers, below).""" wv = wetness(rain, dry, soak) if wet is None else wet ea = soil * (1 - wv) + under(COLUMN if column is None else column) * (1 - (SOAK if soak is None else soak)) return 1 - ea, ea sea_depth = (sea_by_depth * np.array([(a_ + b_) / 2 for a_, b_ in zip(edges[:-1], edges[1:])])[:, None, None]).sum(0) / np.where(1 - near > 0, 1 - near, 1) SREF = gm(sunY[NOW[0] - Y0:NOW[1] - Y0 + 1].mean((0, 1))) # the planet's mean sunlight at the ground, W m-2 def den(sun_sum, months, by='amount'): """What the light's three fluxes, summed over `months` months, are divided by: the planet's mean sunlight (so that a bright place has more light than a dim one), or, for comparison, the place's own sunlight.""" return np.full(np.shape(sun_sum), SREF * months) if by == 'amount' else sun_sum def period(y0, y1, by='amount'): a, b = y0 - Y0, y1 - Y0 + 1 return L3Y[a:b].sum((0, 1)), den(sunY[a:b].sum((0, 1)), 12 * (b - a), by), rainY[a:b].sum((0, 1)), dryY[a:b].sum((0, 1)) def shares(y0, y1, G=None, by='amount', soak=None, column=None): l3, s, rn, dr = period(y0, y1, by) return join(l3, s, *ground(rn, dr, soak, column), G=G) def pool(l3, s, wa, ea, ww, G=None): """A region: its light pooled by sunlight, its ground by area, then joined.""" return join((l3 * ww[..., None]).sum((-3, -2)), (s * ww).sum((-2, -1)), (wa * ww).sum((-2, -1)) / ww.sum(), (ea * ww).sum((-2, -1)) / ww.sum(), G=G) def pooled(mask, y0, y1, G=None, by='amount', soak=None, column=None): l3, s, rn, dr = period(y0, y1, by) return pool(l3, s, *ground(rn, dr, soak, column), W * mask, G=G) ALL = np.ones((NLAT, NLON), bool) def gshares(y0, y1, G=None, by='amount', soak=None, column=None): return pooled(ALL, y0, y1, G, by, soak, column) def placewise(Sbox, mask): """The doshas of a region: the weave made at each box, the boxes averaged by area.""" ww = W * mask return (Sbox * ww[..., None]).sum((0, 1)) / ww.sum() M = {'source': {'product': 'ERA-Interim (ECMWF), as prepared for obs4MIPs by PCMDI', 'via': 'Climate-REF/ref-sample-data, commit 08a1e3a (15 Dec 2025)', 'grid': {'nlat': NLAT, 'nlon': NLON, 'points': NLAT * NLON, 'distinct': (NLAT - 2) * NLON + 2, 'step': 10}, 'years': [Y0, Y1], 'nyears': NY, 'now': list(NOW), 'then': list(THEN)}, 'weave': {k: r(v, 4) for k, v in WEAVE.items()}} # ---------- the plain budget, as a check against published means ---------- a, b = NOW[0] - Y0, NOW[1] - Y0 + 1 rawN = {k: v.reshape(NY, 12, NLAT, NLON)[a:b].mean((0, 1)) for k, v in P['raw'].items()} raw = {k: gm(v) for k, v in rawN.items()} M['budget'] = {k: r(v, 1) for k, v in raw.items()} M['budget']['absorbed'] = r(raw['rsds'] - raw['rsus'], 1) outgo = (pos(P['raw']['lw']) + pos(P['raw']['lh']) + pos(P['raw']['sh'])).reshape(NY, 12, NLAT, NLON)[a:b].mean((0, 1)) paid = (P['air'] + P['fire'] + P['evap']).reshape(NY, 12, NLAT, NLON)[a:b].mean((0, 1)) M['budget']['paidBySun'] = pct(gm(paid) / gm(outgo), 1) # ---------- land, sea and coast ---------- land, sea = near >= 0.999, near <= 0.001 coast = ~land & ~sea M['land'] = {'landShare': pct(gm(land), 1), 'seaShare': pct(gm(sea), 1), 'coastShare': pct(gm(coast), 1), 'boxShare': pct(gm(landfrac), 1), 'landCount': int(land[1:-1].sum()), 'seaCount': int(sea[1:-1].sum()), 'coastCount': int(coast[1:-1].sum())} # for the comparison with Wild et al. every point is counted, as land where most of what it was read from is land for nm, mask in (('budgetLand', near >= 0.5), ('budgetSea', near < 0.5)): ww = W * mask / (W * mask).sum() rr = {k: float((v * ww).sum()) for k, v in rawN.items()} M[nm] = {k: r(v, 0) for k, v in rr.items()} M[nm]['absorbed'] = r(rr['rsds'] - rr['rsus'], 0) # ---------- the ground, today ---------- l3N, sN, rainN, dryN = period(*NOW, by='share') waN, eaN = ground(rainN, dryN) wetN = wetness(rainN, dryN) mm = 86400 * 365.25 / (12 * (NOW[1] - NOW[0] + 1)) # summed monthly rates -> mm a year M['ground'] = {'land': pct(gm(near), 1), 'ice': pct(gm(ice), 1), 'lake': pct(gm(lake), 2), 'soil': pct(gm(soil), 1), 'water': pct(gm(waN), 1), 'earth': pct(gm(eaN), 1), 'wetLand': pct((wetN * soil * W).sum() / (soil * W).sum(), 0), # how much of the world's soil is water, by this rule 'rainLand': r((rainN * mm * soil * W).sum() / (soil * W).sum(), 0), 'dryLand': r((dryN * mm * soil * W).sum() / (soil * W).sum(), 0), 'weight': GROUND, 'column': COLUMN, 'iceThickness': ICE, 'seaDepth': r((sea_depth * (1 - near) * W).sum() / ((1 - near) * W).sum(), 0), # the mean depth of the sea as these boxes read it, m 'seaFloor': pct((under(COLUMN) * (near <= 0.001) * W).sum() / ((near <= 0.001) * W).sum(), 0), # how much of the open sea's column is floor 'seaEarth': pct((under(COLUMN) * (1 - SOAK) * (near <= 0.001) * W).sum() / ((near <= 0.001) * W).sum(), 0), # and how much of its ground is therefore earth 'seaNoFloor': pct(((under(COLUMN) <= 0) * (near <= 0.001) * W).sum() / ((near <= 0.001) * W).sum(), 0), # how much of the open sea is deeper than the column throughout 'seaBelowColumn': pct((sea_by_depth * np.array([float(np.mean(np.linspace(a_, b_, 201) > COLUMN)) for a_, b_ in zip(edges[:-1], edges[1:])])[:, None, None] * sea * W).sum() / (sea_by_depth * sea * W).sum(), 0), # of the open sea's area, how much has its floor below the column (depth taken as even across each class) 'iceFloor': pct(max(0.0, 1 - ICE / COLUMN), 0), 'iceEarth': pct(max(0.0, 1 - ICE / COLUMN) * (1 - SOAK), 0), 'seaByDepth': pct([(sea_by_depth[k] * W).sum() / ((1 - soil - ice) * W).sum() for k in range(len(edges) - 1)], 1), 'depthEdges': edges} # where the counted water and earth are: shares of the planet's ground, by kind col_w = np.array([float(np.mean(np.minimum(1.0, np.linspace(a_, b_, 201) / COLUMN))) for a_, b_ in zip(edges[:-1], edges[1:])]) # how much of the column is sea water, by depth class kw = {'sea': gm((sea_by_depth * col_w[:, None, None]).sum(0)), 'ice': gm(ice * min(1.0, ICE / COLUMN)), 'underWater': gm(under(COLUMN) * SOAK), 'soil': gm(soil * wetN)} ke = {'land': gm(soil * (1 - wetN)), 'seaFloor': gm((sea_by_depth * floor_share(COLUMN)[:, None, None]).sum(0) * (1 - SOAK)), 'underIce': gm(ice * max(0.0, 1 - ICE / COLUMN) * (1 - SOAK))} M['ground']['waterKinds'] = {k: pct(v / sum(kw.values()), 0) for k, v in kw.items()} M['ground']['earthKinds'] = {k: pct(v / sum(ke.values()), 0) for k, v in ke.items()} assert abs(sum(kw.values()) - gm(waN)) < 1e-6 and abs(sum(ke.values()) - gm(eaN)) < 1e-6 lN = l3N / sN[..., None] gl = (l3N * W[..., None]).sum((0, 1)) / (sN * W).sum() M['light'] = {'shares': pct(gl, 1), 'total': pct(gl.sum(), 1), 'intoGround': pct(1 - gl.sum(), 1), 'reference': r(SREF, 1), 'air': {'moving': MOVING, 'meanWind': r(WBAR, 1), 'meanHeat': r(ABAR, 1), # how far each half strays from place to place (standard deviation over the surface of its ten-year mean, W m-2), over the planet and over land 'spreadHeat': r(np.sqrt(gm((HEAT[n0_:n1_].mean(0) - ABAR) ** 2)), 1), 'spreadMove': r(np.sqrt(gm((MOVE[n0_:n1_].mean(0) - ABAR) ** 2)), 1), 'windLand': r((windM[n0_:n1_].mean(0) * W * (near >= 0.999)).sum() / (W * (near >= 0.999)).sum(), 1), 'windSea': r((windM[n0_:n1_].mean(0) * W * (near <= 0.001)).sum() / (W * (near <= 0.001)).sum(), 1)}} amount = l3N.sum(-1) / (SREF * 12 * (NOW[1] - NOW[0] + 1)) # the light of each box, in hundreds: 1 would be all of the planet's mean sunlight kept M['light']['amount'] = {'min': pct(amount.min(), 0), 'max': pct(amount.max(), 0)} # ---------- today's map ---------- E = shares(*NOW) # lat, lon, 5 gE = gshares(*NOW) S = strengths(E) # the weave at each box D = norm(S) gS = strengths(gE) # the planet taken whole gD = norm(gS) ratio = S / gS # each dosha against the strength it has on the planet as a whole lead = S.argmax(-1) M['now'] = {'elements': pct(gE, 1), 'elements2': pct(gE, 2), 'strengths': r(gS, 3), 'doshas': pct(gD, 1), 'doshasByPlace': pct(norm(placewise(S, ALL)), 1), 'sums': pct([gE[0] + gE[1], gE[2], gE[3] + gE[4]], 1)} icy = ice >= 0.999 M['now']['ice'] = {'area': pct(W[icy].sum(), 1), 'elements': pct(pooled(icy, *NOW), 0), 'doshas': pct(norm(placewise(S, icy)), 0), 'lead': pct([(W * icy)[lead == k].sum() / (W * icy).sum() for k in range(3)], 0)} mixed = (near > 0.001) & (near < 0.999) M['now']['coast'] = {'area': pct(W[mixed].sum(), 1), 'elements': pct(pooled(mixed, *NOW), 0), 'doshas': pct(norm(placewise(S, mixed)), 0), 'lead': pct([(W * mixed)[lead == k].sum() / (W * mixed).sum() for k in range(3)], 0)} for nm, mask in (('land', land), ('sea', sea)): M['now'][nm] = {'elements': pct(pooled(mask, *NOW), 1), 'doshas': pct(norm(placewise(S, mask)), 1), 'area': pct((W * mask).sum(), 1)} def lead_area(mask): ww = W * mask return pct([ww[lead == k].sum() / ww.sum() for k in range(3)], 1) M['now']['lead'] = {'area': lead_area(ALL), 'landArea': lead_area(land), 'seaArea': lead_area(sea)} # by latitude band, and by longitude band (30 degrees wide, 60 S to 60 N) M['now']['bands'] = [] for j in range(NLAT - 1, -1, -1): mask = np.zeros((NLAT, NLON), bool) mask[j] = True M['now']['bands'].append({'lat': int(lat[j]), 'elements': pct(pooled(mask, *NOW), 1), 'sun': r(sunY[a:b, :, j].mean(), 0), 'land': pct(near[j].mean(), 0), 'doshas': pct(norm(placewise(S, mask)), 1)}) M['now']['lonBands'] = [] jj = (lat >= -60) & (lat <= 60) for i0 in range(0, NLON, 3): mask = np.zeros((NLAT, NLON), bool) mask[np.ix_(jj, [(i0 - 1) % NLON, i0, (i0 + 1) % NLON])] = True M['now']['lonBands'].append({'lon': int(lon[i0]), 'elements': pct(pooled(mask, *NOW), 1), 'doshas': pct(norm(placewise(S, mask)), 1), 'land': pct((landfrac * W * mask).sum() / (W * mask).sum(), 0)}) # named places (the box whose centre is nearest) def kind(j, i): return 'land' if land[j, i] else 'sea' if sea[j, i] else 'coast' def place(name, la, lo): j, i = int(round((la + 90) / 10)), int(round((lo % 360) / 10)) % NLON return {'name': name, 'lat': int(lat[j]), 'lon': int(lon[i]) if lon[i] <= 180 else int(lon[i] - 360), 'kind': kind(j, i), 'landShare': pct(near[j, i], 0), 'wet': pct(wetN[j, i], 0), 'light': pct(lN[j, i], 0), 'lightAmount': pct(amount[j, i], 0), 'sun': r(sN[j, i] / (12 * (NOW[1] - NOW[0] + 1)), 0), 'elements': pct(E[j, i], 1), 'doshas': pct(D[j, i], 0), 'lead': 'VPK'[lead[j, i]]} M['now']['places'] = [place(*p) for p in [ ('Sahara', 20, 10), ('Arabia', 20, 50), ('Central Australia', -20, 130), ('Kalahari', -20, 20), ('Central India', 20, 80), ('Amazon', 0, -60), ('Congo', 0, 20), ('Borneo', 0, 110), ('North Atlantic', 50, -30), ('North Pacific', 40, -170), ('Equatorial Pacific', 0, -120), ('Southern Ocean', -50, 90), ('Trade-wind Atlantic', 20, -40), ('Siberia', 60, 100), ('Central Europe', 50, 10), ('Great Plains', 40, -100), ('Greenland', 70, -40), ('Antarctica', -80, 60), ('Arctic Ocean', 80, 180)]] # ---------- the three lands, told apart by measured rain (GPCP) ---------- rain = load(data_dir, GPCP, 'pr', f'{NOW[0]}-01', f'{NOW[1]}-12')[0].mean(0) * 86400 * 365.25 # mm a year temp = load(data_dir, ERA, 'ts', f'{NOW[0]}-01', f'{NOW[1]}-12')[0].mean(0) - 273.15 wind = load(data_dir, ERA, 'sfcWind', f'{NOW[0]}-01', f'{NOW[1]}-12')[0].mean(0) temperate = (np.abs(lat)[:, None] <= 60) & ALL def lands(lo, hi): row = {'dryBelow': lo, 'wetAbove': hi} for nm, mask in (('dry', land & temperate & (rain < lo)), ('mid', land & temperate & (rain >= lo) & (rain <= hi)), ('wet', land & temperate & (rain > hi))): ww = W * mask dd = pct(norm(placewise(S, mask)), 0) box = D[mask] row[nm] = {'count': int(mask.sum()), 'rain': r((rain * ww).sum() / ww.sum(), 0), 'wet': pct((wetN * ww).sum() / ww.sum(), 0), 'elements': pct(pooled(mask, *NOW), 0), 'doshas': dd, 'spread': r(max(dd) - min(dd), 0), 'lead': pct([ww[lead == k].sum() / ww.sum() for k in range(3)], 0), 'median': pct(np.median(box, 0), 0), 'pooledDoshas': pct(norm(strengths(pooled(mask, *NOW))), 0)} return row M['now']['lands'] = lands(250, 1500) M['now']['landsAlt'] = [lands(150, 1000), lands(400, 2000)] # the three doshas against rain, temperature and wind, box by box def cc(x, mask): return [r(corr(D[..., k][mask], x[mask], W[mask]), 2) for k in range(3)] lt, st = land & temperate, sea & temperate M['now']['against'] = {'landRain': cc(np.log(rain + 1), lt), 'landTemp': cc(temp, lt), 'landWind': cc(wind, lt), 'seaRain': cc(np.log(rain + 1), st), 'seaTemp': cc(temp, st), 'seaWind': cc(wind, st), 'landCount': int(lt.sum()), 'seaCount': int(st.sum()), 'windLand': r((wind * W * land).sum() / (W * land).sum(), 1), 'windSea': r((wind * W * sea).sum() / (W * sea).sum(), 1)} # ---------- CERES as a check on the radiation ---------- C = {v: load(data_dir, CERES, v, '2001-01', '2018-12')[0] for v in CERES_VARS} Fc = {v: load(data_dir, ERA, v, '2001-01', '2018-12')[0] for v in ERA_VARS} chk = {} for nm, fe, fc in (('rsds', Fc['rsds'], C['rsds']), ('rsus', Fc['rsus'], C['rsus']), ('lwnet', Fc['rlus'] - Fc['rlds'], C['rlus'] - C['rlds'])): chk[nm] = {'era': r(gm(fe.mean(0)), 1), 'ceres': r(gm(fc.mean(0)), 1), 'pattern': r(corr(fe.mean(0), fc.mean(0), W), 3)} Fmix = dict(Fc) Fmix.update(C) Pm, Pe = split(Fmix), split(Fc) rainC = pos(load(data_dir, GPCP, 'pr', '2001-01', '2018-12')[0]).sum(0) def five_from(Pp, G=None): l3m = light3(Pp) l3m[..., 1] = air_of(l3m[..., 1], MOVE[(2001 - Y0) * 12:(2018 - Y0 + 1) * 12]) l3, s = l3m.sum(0), den(Pp['rsds'].sum(0), Pp['rsds'].shape[0]) wa, ea = ground(rainC, Pp['dry'].sum(0)) return join(l3, s, wa, ea, G=G), pool(l3, s, wa, ea, W, G=G) Em, gEm = five_from(Pm) Ee, gEe = five_from(Pe) chk.update({'sharesEra': pct(gEe, 1), 'sharesCeres': pct(gEm, 1), 'doshasEra': pct(norm(strengths(gEe)), 1), 'doshasCeres': pct(norm(strengths(gEm)), 1), 'sameLead': pct(W[strengths(Em).argmax(-1) == strengths(Ee).argmax(-1)].sum(), 0), 'years': [2001, 2018]}) M['ceres'] = chk # ---------- how much the reading leans on the proportions, and on the source of the radiation ---------- M['now']['recipes'] = {} for nm, Wv in (('caraka', WEAVE), ('tastes', TASTES), ('equal', EQUAL)): lv = strengths(E, Wv).argmax(-1) M['now']['recipes'][nm] = {'doshas': pct(norm(strengths(gE, Wv)), 1), 'doshasByPlace': pct(norm(placewise(strengths(E, Wv), ALL)), 1), 'doshasSatellite': pct(norm(strengths(gEm, Wv)), 1), 'lead': pct([W[lv == k].sum() for k in range(3)], 1)} # ---------- how much it leans on the weight given to the ground against the light ---------- M['now']['weights'] = {} for nm, G, by, sk, col in (('half', 0.5, 'amount', None, None), ('one', 1.0, 'amount', None, None), ('double', 2.0, 'amount', None, None), ('fixed', 'fixed', 'amount', None, None), ('byShare', 1.0, 'share', None, None), ('soakThird', 1.0, 'amount', 1 / 3, None), ('soakTwoThirds', 1.0, 'amount', 2 / 3, None), ('surface', 1.0, 'amount', None, 0.0), ('column200', 1.0, 'amount', None, 200.0), ('column1000', 1.0, 'amount', None, 1000.0), ('column2000', 1.0, 'amount', None, 2000.0), ('column6000', 1.0, 'amount', None, 6000.0)): Eg, gg = shares(*NOW, G=G, by=by, soak=sk, column=col), gshares(*NOW, G=G, by=by, soak=sk, column=col) Sg = strengths(Eg) lg = Sg.argmax(-1) lt_, wt_ = land & (np.abs(lat)[:, None] <= 60) & (rain < 250), land & (np.abs(lat)[:, None] <= 60) & (rain > 1500) M['now']['weights'][nm] = {'elements': pct(gg, 1), 'doshas': pct(norm(strengths(gg)), 1), 'doshasByPlace': pct(norm(placewise(Sg, ALL)), 1), 'lead': pct([W[lg == k].sum() for k in range(3)], 1), 'leadSea': pct([(W * sea)[lg == k].sum() / (W * sea).sum() for k in range(3)], 0), 'sea': pct(norm(placewise(Sg, sea)), 0), 'ice': pct(norm(placewise(Sg, ice >= 0.999)), 0), 'seaFloor': pct((under(COLUMN if col is None else col) * sea * W).sum() / (sea * W).sum(), 0), 'seaEarth': pct(((Eg[..., 4] / (Eg[..., 3] + Eg[..., 4])) * sea * W).sum() / (sea * W).sum(), 0), 'seaNoFloor': pct(((under(COLUMN if col is None else col) <= 0) * sea * W).sum() / (sea * W).sum(), 0), 'leadLand': pct([(W * land)[lg == k].sum() / (W * land).sum() for k in range(3)], 0), 'dry': pct(norm(placewise(Sg, lt_)), 0), 'mid': pct(norm(placewise(Sg, land & (np.abs(lat)[:, None] <= 60) & ~lt_ & ~wt_)), 0), 'wet': pct(norm(placewise(Sg, wt_)), 0)} # ---------- the step of the accounting: the light split made on whole years, not months ---------- Fyear = {v: F[v].reshape(NY, 12, NLAT, NLON)[a:b].mean((0, 1)) for v in ERA_VARS} Py = split(Fyear) l3y_ = light3(Py) warm_y = pct((l3y_[..., 1] * W).sum() / (Py['rsds'] * W).sum(), 1) # the warming of the air alone, by the yearly accounting l3y_[..., 1] = air_of(l3y_[..., 1], MOVE[n0_:n1_].mean(0)) M['now']['yearlyAccounting'] = {'elements': pct(pool(l3y_, den(Py['rsds'], 1), waN, eaN, W), 1), 'light': pct((l3y_ * W[..., None]).sum((0, 1)) / (Py['rsds'] * W).sum(), 1), 'warming': warm_y} # ---------- the year ---------- cm = L3Y[a:b].sum(0) # 12, lat, lon, 3 cs = sunY[a:b].sum(0) # 12, lat, lon dark = cs / (b - a) <= DARK # polar night, near enough; in such a month a place is its ground alone wa12, ea12 = ground(rainY[a:b].sum(0), dryY[a:b].sum(0)) # the ground month by month: wetter when rain outruns the drying wet12 = wetness(rainY[a:b].sum(0), dryY[a:b].sum(0)) Em12 = join(cm, den(cs, b - a), wa12, ea12) Sm12 = strengths(Em12) # 12, lat, lon, 3: the weave at each box in each month M['year'] = {'darkMonths': pct((dark * W[None]).sum() / 12, 1)} def season_string(rel): return ''.join('VPK'[k] if rel[m, k] > 1 else '-' for m, k in enumerate(rel.argmax(-1))) def cyc_order(peaks): v, p, k = peaks return 'VKP' if (k - v) % 12 < (p - v) % 12 else 'VPK' def calendar(mask, name): ww = W * mask s12 = (Sm12 * ww[None, ..., None]).sum((1, 2)) / ww.sum() # region: boxes averaged own = (S * ww[..., None]).sum((0, 1)) / ww.sum() rel = s12 / np.where(own > 0, own, np.nan) # a dosha with no measure at all in a region (kapha over open sea) has no calendar peaks = [int(rel[:, k].argmax()) for k in range(3)] out = {'name': name, 'count': int(mask.sum()), 'elements': [pct(x, 0) for x in pool(cm, den(cs, b - a), wa12, ea12, ww)], 'sun': r((cs * ww[None]).sum((1, 2)) / ww.sum() / (b - a), 0), 'wet': pct((wet12 * (soil * ww)[None]).sum((1, 2)) / max((soil * ww).sum(), 1e-12), 0), 'ratio': [r(x, 2) for x in rel], 'season': season_string(rel), 'above': [''.join('1' if rel[m, k] > 1 else '0' for m in range(12)) for k in range(3)], 'has': [bool(own[k] > 0) for k in range(3)], 'peak': peaks, 'peakRatio': [r(rel[peaks[k], k], 2) for k in range(3)], 'low': [r(rel[:, k].min(), 2) for k in range(3)], 'order': cyc_order(peaks)} # box by box if mask.sum() > 1: relb = Sm12[:, mask] / np.where(S[mask] > 0, S[mask], 1)[None] # 12, boxes, 3 wb = ww[mask] pk = relb.argmax(0) # boxes, 3 orders = np.array([cyc_order(x) for x in pk]) out['boxes'] = {'order': {o: pct(wb[orders == o].sum() / wb.sum(), 0) for o in ('VKP', 'VPK')}, 'orderCount': {o: int((orders == o).sum()) for o in ('VKP', 'VPK')}, 'peakMonthMedian': [int(np.median(pk[:, k])) for k in range(3)], # in how many boxes vata's highest month falls before, in, and after kapha's (plain months of the calendar year) 'vataVsKapha': [int((pk[:, 0] < pk[:, 2]).sum()), int((pk[:, 0] == pk[:, 2]).sum()), int((pk[:, 0] > pk[:, 2]).sum())]} return out band = lambda lo, hi: ((lat >= lo) & (lat <= hi))[:, None] & np.ones((1, NLON), bool) one = lambda la, lo: (lat[:, None] == la) & (lon[None, :] == lo % 360) M['year']['calendars'] = { 'nLand': calendar(band(30, 50) & land, 'land, 30 to 50 N'), 'nSea': calendar(band(30, 50) & sea, 'sea, 30 to 50 N'), 'india': calendar(one(20, 80), 'central India, 20 N 80 E'), 'sahel': calendar((lat[:, None] == 10) & ((lon[None, :] <= 30) | (lon[None, :] >= 350)) & land, 'the Sahel, 10 N, 10 W to 30 E'), 'plains': calendar(one(40, -100), 'Great Plains, 40 N 100 W'), 'europe': calendar(one(50, 10), 'central Europe, 50 N 10 E'), 'sLand': calendar(band(-40, -20) & land, 'land, 20 to 40 S'), 'tropLand': calendar(band(-10, 10) & land, 'land, 10 S to 10 N'), } # how much of the surface keeps one leading dosha through the year leadm = np.where(Sm12.max(-1) > 0, Sm12.argmax(-1), -1) # -1 in a month with nothing woven at all (open sea in polar night) nlead = np.array([[len(set(leadm[:, j, i]) - {-1}) for i in range(NLON)] for j in range(NLAT)]) for nm, mask in (('', ALL), ('Land', land), ('Sea', sea)): M['year']['oneLeadAllYear' + nm] = pct((W * mask)[nlead <= 1].sum() / (W * mask).sum(), 0) pitta_all = np.array([[set(leadm[:, j, i]) - {-1} == {1} for i in range(NLON)] for j in range(NLAT)]) M['year']['pittaAllYearSea'] = pct((W * sea)[pitta_all].sum() / (W * sea).sum(), 0) # of the open sea, where pitta leads in every month # how far a month stands from the same place's year dm = np.abs(Em12 - E[None]).sum(-1) / 2 M['year']['monthFromYear'] = pct((dm * W[None]).sum() / 12, 1) # ---------- two tests: a fire that knows how warm the place is, and a soil that remembers ---------- # Neither is used anywhere else. They are run beside the main count to see what each would change. tsY = load(data_dir, ERA, 'ts', f'{Y0}-01', f'{Y1}-12')[0].reshape(NY, 12, NLAT, NLON) # surface temperature, K TBAR = gm(tsY[a:b].mean((0, 1))) # the planet's mean surface temperature today sden = np.full((NLAT, NLON), SREF * 12 * (b - a)) rain12, dry12 = rainY[a:b].sum(0), dryY[a:b].sum(0) def cal_core(mask, Sm, Sown): ww = W * mask s12 = (Sm * ww[None, ..., None]).sum((1, 2)) / ww.sum() own = (Sown * ww[..., None]).sum((0, 1)) / ww.sum() rel = s12 / np.where(own > 0, own, np.nan) peaks = [int(np.nanargmax(rel[:, k])) if own[k] > 0 else 0 for k in range(3)] return {'count': int(mask.sum()), 'ratio': [[None if np.isnan(v) else r(v, 2) for v in x] for x in rel], 'above': [''.join('1' if rel[m, k] > 1 else '0' for m in range(12)) for k in range(3)], 'has': [bool(own[k] > 0) for k in range(3)], 'peak': peaks, 'peakRatio': [r(rel[peaks[k], k], 2) if own[k] > 0 else None for k in range(3)], 'low': [r(np.nanmin(rel[:, k]), 2) if own[k] > 0 else None for k in range(3)]} CAL = {'nLand': band(30, 50) & land, 'india': one(20, 80), 'sahel': (lat[:, None] == 10) & ((lon[None, :] <= 30) | (lon[None, :] >= 350)) & land, 'sLand': band(-40, -20) & land, 'plains': one(40, -100), 'europe': one(50, 10)} def test(l3m, wet12v=None, wetNv=None): """The main count with the monthly light, or the wetness of the soil, replaced.""" l3y = l3m.reshape(NY, 12, NLAT, NLON, 3) l3n = l3y[a:b].sum((0, 1)) waV, eaV = ground(rainN, dryN, wet=wetNv) Ev, gEv = join(l3n, sden, waV, eaV), pool(l3n, sden, waV, eaV, W) Sv = strengths(Ev) Dv, lv = norm(Sv), Sv.argmax(-1) S0 = strengths(join(l3n, sden, *ground(rainN, dryN, column=0.0, wet=wetNv))) # the same with the ground counted at its surface only ice0 = placewise(S0, icy) iceMain = placewise(strengths(shares(*NOW, column=0.0)), icy) Smv = strengths(join(l3y[a:b].sum(0), den(cs, b - a), *ground(rain12, dry12, wet=wet12v))) mid_ = lt & (rain >= 250) & (rain <= 1500) out = {'elements': pct(gEv, 1), 'doshas': pct(norm(strengths(gEv)), 1), 'doshasByPlace': pct(norm(placewise(Sv, ALL)), 1), 'lead': pct([W[lv == k].sum() for k in range(3)], 1), 'leadLand': pct([(W * land)[lv == k].sum() / (W * land).sum() for k in range(3)], 0), 'sea': pct(norm(placewise(Sv, sea)), 0), 'ice': pct(norm(placewise(Sv, icy)), 0), 'leadSea': pct([(W * sea)[lv == k].sum() / (W * sea).sum() for k in range(3)], 0), 'seaSouth': pct(norm(placewise(Sv, sea & band(-60, -50))), 0), 'seaTropics': pct(norm(placewise(Sv, sea & band(-20, 20))), 0), # the open sea at 50 and 60 south, and within 20 degrees of the equator 'iceSurface': pct(norm(ice0), 0), 'iceSurfaceWoven': pct(ice0.sum() / iceMain.sum(), 0), # the ice with the ground counted at its surface, and how much is woven there at all against the main count 'dry': pct(norm(placewise(Sv, lt & (rain < 250))), 0), 'mid': pct(norm(placewise(Sv, mid_)), 0), 'wet': pct(norm(placewise(Sv, lt & (rain > 1500))), 0), 'dryLead': pct([(W * (lt & (rain < 250)))[lv == k].sum() / (W * (lt & (rain < 250))).sum() for k in range(3)], 0), 'landTemp': [r(corr(Dv[..., k][lt], temp[lt], W[lt]), 2) for k in range(3)], 'landRain': [r(corr(Dv[..., k][lt], np.log(rain + 1)[lt], W[lt]), 2) for k in range(3)], 'landWind': [r(corr(Dv[..., k][lt], wind[lt], W[lt]), 2) for k in range(3)], 'seaTemp': [r(corr(Dv[..., k][st], temp[st], W[st]), 2) for k in range(3)], 'seaWind': [r(corr(Dv[..., k][st], wind[st], W[st]), 2) for k in range(3)], 'calendars': {k: cal_core(m_, Smv, Sv) for k, m_ in CAL.items()}, 'places': {nm: pct(Dv[int(round((la + 90) / 10)), int(round((lo % 360) / 10)) % NLON], 0) for nm, la, lo in (('Sahara', 20, 10), ('Central India', 20, 80), ('Amazon', 0, -60), ('Central Europe', 50, 10), ('Great Plains', 40, -100), ('Siberia', 60, 100), ('Greenland', 70, -40), ('Antarctica', -80, 60))}} return out, Ev, gEv def with_fire(h): l3 = L3.copy() l3[..., 2] = l3[..., 2] * h.reshape(NY * 12, NLAT, NLON) return l3 STORE = {'capacity': 150, 'fullPowerAbove': 0.75, 'stepsPerMonth': 30, 'sizesTried': [75, 150, 300]} def bucket(cap, sub=STORE['stepsPerMonth']): """A soil that holds up to `cap` millimetres and carries what it holds from month to month. Rain fills it; the sun dries it at its full power while the soil is more than three-quarters full and in proportion below that (Manabe 1969); what overflows runs off. Returns the wetness of the soil month by month: half of how full it is.""" mmm = 86400 * 30.4375 Pm, Dm = rainM * mmm, P['dry'] * mmm Sb, out = np.full((NLAT, NLON), cap / 2.0), np.zeros((NY * 12, NLAT, NLON)) for _ in (0, 1): # the first pass through the record only sets the starting state for m in range(NY * 12): acc = 0.0 for _s in range(sub): Sb = np.clip(Sb + Pm[m] / sub - Dm[m] / sub * np.minimum(1.0, Sb / (STORE['fullPowerAbove'] * cap)), 0, cap) acc = acc + Sb out[m] = acc / sub / cap return SOAK * out.reshape(NY, 12, NLAT, NLON) warm = np.clip(tsY - 273.15, 0, None) / (TBAR - 273.15) # nothing at freezing, one at the planet's mean emit = (tsY / TBAR) ** 4 # what a surface at that temperature radiates, against the planet's mean T = {} T['main'], _, _ = test(L3) T['warm'], Ewarm, gEwarm = test(with_fire(warm)) T['emit'], _, _ = test(with_fire(emit)) wetB = {cap: bucket(cap)[a:b].mean(0) for cap in (75, 150, 300)} # 12, lat, lon for cap in (75, 150, 300): T['bucket%d' % cap], Eb, gEb = test(L3, wetB[cap], wetB[cap].mean(0)) T['bucket%d' % cap]['wetLand'] = pct((wetB[cap].mean(0) * soil * W).sum() / (soil * W).sum(), 0) T['bucket%d' % cap]['wetAgrees'] = r(corr(wetB[cap].mean(0)[lt], wetN[lt], W[lt]), 2) T['bucket%d' % cap]['wetByMonth'] = {k: pct((wetB[cap] * (soil * W * m_)[None]).sum((1, 2)) / (soil * W * m_).sum(), 0) for k, m_ in CAL.items()} T['both'], _, _ = test(with_fire(warm), wetB[150], wetB[150].mean(0)) # the air counted by the heat handed to it alone, and by its movement alone, and at steps between def with_moving(m): l3 = L3.copy() l3[..., 1] = air_of(HEAT, MOVE, m) return l3 T['heatOnly'], Eheat, gEheat = test(with_moving(0.0)) T['moveOnly'], Emove, gEmove = test(with_moving(1.0)) T['moving'] = {} for m_ in (0.25, 0.75, 0.9, 0.95): t_, _, _ = test(with_moving(m_)) T['moving']['%g' % m_] = {'dry': t_['dry'], 'dryLead': t_['dryLead'], 'sea': t_['sea'], 'india': t_['calendars']['india']['above'][0], 'indiaPeak': t_['calendars']['india']['peak'][0], 'indiaJulyOverMay': r(t_['calendars']['india']['ratio'][6][0] / t_['calendars']['india']['ratio'][4][0], 2)} T['wind'] = {k: [r(float((windM.reshape(NY, 12, NLAT, NLON)[a:b].mean(0)[m] * W * m_).sum() / (W * m_).sum()), 1) for m in range(12)] for k, m_ in CAL.items()} T['heat'] = {k: [r(float((HEAT.reshape(NY, 12, NLAT, NLON)[a:b].mean(0)[m] * W * m_).sum() / (W * m_).sum()), 0) for m in range(12)] for k, m_ in CAL.items()} # the main count's own monthly wetness, with the year's wetness taken as the mean of its months (as the store's is), to tell what the store adds from what the averaging does T['monthMean'], _, _ = test(L3, wet12, wet12.mean(0)) T['storeOnly'], _, _ = test(L3, wetB[150], None) # the store's months held against the year as the main count reckons it: memory with the measure left alone T['warm']['iceWarmestMonth'] = r(float((tsY[a:b] - 273.15)[:, :, icy].max()), 1) # the warmest month on any box that is wholly ice, deg C T['warm']['fire'] = {'frozenMonths': pct((W[None, None] * (tsY[a:b] <= 273.15)).sum() / (12 * (b - a)), 0), 'hottest': r(float(warm[a:b].mean((0, 1)).max()), 1), 'sahara': r(float(warm[a:b].mean((0, 1))[int(round((20 + 90) / 10)), 1]), 1)} T['meanTemp'] = r(TBAR - 273.15, 1) T['store'] = STORE T['warmth'] = {'frozenArea': pct((W * (tsY[a:b].mean((0, 1)) <= 273.15)).sum(), 0), 'max': r(float(warm[a:b].mean((0, 1)).max()), 2)} M['tests'] = T # ---------- forty years ---------- waY, eaY = ground(rainY.sum(1), dryY.sum(1)) EY = join(L3Y.sum(1), den(sunY.sum(1), 12), waY, eaY) # year, lat, lon, 5 gY = np.array([gshares(y, y) for y in range(Y0, Y1 + 1)]) SY = strengths(EY) gDY = norm(strengths(gY)) yrs = np.arange(NY) def trend(x): c = np.polyfit(yrs, x, 1) resid = x - np.polyval(c, yrs) se = np.sqrt((resid ** 2).sum() / (NY - 2) / ((yrs - yrs.mean()) ** 2).sum()) return c[0] * 10, se * 10, resid.std() eT, eN = gshares(*THEN), gshares(*NOW) M['forty'] = {'elementsThen': pct(eT, 3), 'elementsNow': pct(eN, 3), 'doshasThen': pct(norm(strengths(eT)), 3), 'doshasNow': pct(norm(strengths(eN)), 3), 'elementTrend': [r(trend(gY[:, k] * 100)[0], 2) for k in range(5)], 'elementTrendSE': [r(trend(gY[:, k] * 100)[1], 2) for k in range(5)], 'elementWobble': [r(trend(gY[:, k] * 100)[2], 2) for k in range(5)], 'doshaTrend': [r(trend(gDY[:, k] * 100)[0], 2) for k in range(3)], 'doshaWobble': [r(trend(gDY[:, k] * 100)[2], 2) for k in range(3)]} # the fire that knows cold (one of the two tests above), through the forty years: does it see the planet warm? l3wY = with_fire(warm).reshape(NY, 12, NLAT, NLON, 3) gYw = np.array([pool(l3wY[y].sum(0), np.full((NLAT, NLON), SREF * 12), *ground(rainY[y].sum(0), dryY[y].sum(0)), W) for y in range(NY)]) tsG = np.array([gm(tsY[y].mean(0)) for y in range(NY)]) - 273.15 windG = np.array([gm(windM.reshape(NY, 12, NLAT, NLON)[y].mean(0)) for y in range(NY)]) M['forty']['wind'] = {'then': r(windG[:10].mean(), 2), 'now': r(windG[-10:].mean(), 2), 'trend': r(trend(windG)[0], 3), 'trendSE': r(trend(windG)[1], 3)} M['forty']['lightTotal'] = {'then': r(gY[:10, :3].sum(1).mean() / (1 - gY[:10, :3].sum(1).mean()) * 100, 2), 'now': r(gY[-10:, :3].sum(1).mean() / (1 - gY[-10:, :3].sum(1).mean()) * 100, 2)} # the light's count against the ground's hundred M['tests']['warm']['forty'] = {'fireThen': r(gYw[:10, 2].mean() * 100, 2), 'fireNow': r(gYw[-10:, 2].mean() * 100, 2), 'fireTrend': r(trend(gYw[:, 2] * 100)[0], 2), 'fireTrendSE': r(trend(gYw[:, 2] * 100)[1], 2), 'fireWobble': r(trend(gYw[:, 2] * 100)[2], 2), 'tempThen': r(tsG[:10].mean(), 2), 'tempNow': r(tsG[-10:].mean(), 2), 'tempTrend': r(trend(tsG)[0], 2), 'tempTrendSE': r(trend(tsG)[1], 2)} rawY = {k: np.array([gm(v.reshape(NY, 12, NLAT, NLON)[y].mean(0)) for y in range(NY)]) for k, v in P['raw'].items()} M['forty']['fluxThen'] = {k: r(v[:10].mean(), 1) for k, v in rawY.items()} M['forty']['fluxNow'] = {k: r(v[-10:].mean(), 1) for k, v in rawY.items()} # the ground's two terms year by year: rain on the world's soil and the sun's power to dry it (mm a year), and the wetness they give sw = soil * W / (soil * W).sum() rainLY, dryLY = (rainY.sum(1) * sw).sum((1, 2)) * 86400 * 365.25 / 12, (dryY.sum(1) * sw).sum((1, 2)) * 86400 * 365.25 / 12 wetLY = (wetness(rainY.sum(1), dryY.sum(1)) * sw).sum((1, 2)) * 100 lightY = np.array([(L3Y[y].sum(0) * W[..., None]).sum((0, 1)) / (sunY[y].sum(0) * W).sum() for y in range(NY)]) * 100 M['forty']['ground'] = {'rainThen': r(rainLY[:10].mean(), 0), 'rainNow': r(rainLY[-10:].mean(), 0), 'dryThen': r(dryLY[:10].mean(), 0), 'dryNow': r(dryLY[-10:].mean(), 0), 'wetThen': r(wetLY[:10].mean(), 1), 'wetNow': r(wetLY[-10:].mean(), 1), 'wetTrend': r(trend(wetLY)[0], 2), 'wetTrendSE': r(trend(wetLY)[1], 2), 'wetWobble': r(trend(wetLY)[2], 2), 'rainTrend': r(trend(rainLY)[0], 1), 'rainTrendSE': r(trend(rainLY)[1], 1), 'dryTrend': r(trend(dryLY)[0], 1), 'dryTrendSE': r(trend(dryLY)[1], 1)} M['forty']['light'] = {'then': r(lightY[:10].mean(0), 2), 'now': r(lightY[-10:].mean(0), 2), 'trend': [r(trend(lightY[:, k])[0], 2) for k in range(3)], 'trendSE': [r(trend(lightY[:, k])[1], 2) for k in range(3)]} # constitution: the leading dosha of each place, year by year Wn = W leadY = SY.argmax(-1) always = (leadY == leadY[0]).all(0) flips = (leadY[1:] != leadY[:-1]) for nm, mask in (('', ALL), ('Land', land), ('Sea', sea)): M['forty']['sameLeadAlways' + nm] = pct((Wn * mask)[always].sum() / (Wn * mask).sum(), 0) M['forty']['flipPerYear'] = pct((flips * Wn).sum() / (NY - 1) / Wn.sum(), 1) srt = np.sort(norm(strengths(shares(Y0, Y1))), -1) margin = srt[..., -1] - srt[..., -2] M['forty']['flipsNearBoundary'] = pct((flips * (Wn * (margin < 0.05))).sum() / (flips * Wn).sum(), 0) M['forty']['areaNearBoundary'] = pct(Wn[margin < 0.05].sum() / Wn.sum(), 0) M['forty']['seaNearBoundary'] = pct((Wn * sea)[margin < 0.05].sum() / Wn[margin < 0.05].sum(), 0) # how much of that boundary country is open sea M['forty']['flipsAtSea'] = pct((flips * (Wn * sea)).sum() / (flips * Wn).sum(), 0) # and how much of the changing is done there M['forty']['thenNowSame'] = pct(Wn[strengths(shares(*THEN)).argmax(-1) == strengths(shares(*NOW)).argmax(-1)].sum() / Wn.sum(), 0) la_y = np.array([[Wn[leadY[y] == k].sum() for k in range(3)] for y in range(NY)]) # shares of the whole surface M['forty']['leadAreaThen'] = pct(la_y[:10].mean(0), 1) M['forty']['leadAreaNow'] = pct(la_y[-10:].mean(0), 1) dist = lambda A, B: float((np.abs(A - B).sum(-1) / 2 * W).sum()) E40 = shares(Y0, Y1) M['forty']['yearFromMean'] = pct(np.mean([dist(EY[y], E40) for y in range(NY)]), 1) M['forty']['yearFromPrev'] = pct(np.mean([dist(EY[y], EY[y - 1]) for y in range(1, NY)]), 1) M['forty']['thenFromNow'] = pct(dist(shares(*THEN), shares(*NOW)), 1) M['forty']['placeFromPlace'] = pct(np.mean([dist(E40, np.roll(E40, s, 1)) for s in range(1, NLON)]), 1) # ---------- the dial: chains of chances counted from the measured years ---------- def chain(states, wts, lag=1, valid=None): Pm = np.zeros((3, 3)) for t in range(states.shape[0] - lag): v = wts > 0 if valid is None else (valid[t] & valid[t + lag] & (wts > 0)) np.add.at(Pm, (states[t][v], states[t + lag][v]), wts[v]) return Pm def rows(Cm): return Cm / Cm.sum(1, keepdims=True) def resting(Pn): vals, vecs = np.linalg.eig(Pn.T) v = np.real(vecs[:, np.argmax(np.real(vals))]) return v / v.sum() def reverse(Pn): # the same chain read backward: Bayes' rule with the chain's own resting frequencies pi = resting(Pn) return (Pn * pi[:, None]).T / pi[:, None] def same_record(states, k, wts, valid=None): Ck = chain(states, wts, k, valid) return float(np.trace(Ck) / Ck.sum()) def same_chain(Pn, freq, k): return float((freq * np.diag(np.linalg.matrix_power(Pn, k))).sum()) def back_record(states, k, wts): return rows(chain(states, wts, k).T) # of the places in each state now, the shares in each state k years before def horizon(Pn, tol=0.01): R, rest = reverse(Pn), resting(Pn) for k in range(1, 400): if np.abs(np.linalg.matrix_power(R, k) - rest[None, :]).max() < tol: return k return None lags = [1, 2, 3, 4, 5, 7, 10, 15, 20, 30] KS = list(range(1, 11)) + [15, 20] def describe(states, wts, digits): C1 = chain(states, wts) Pn, freq = rows(C1), C1.sum(1) / C1.sum() ev = np.sort(np.abs(np.linalg.eigvals(Pn)))[::-1] return Pn, {'P': [r(x, digits) for x in Pn], 'frequencies': pct(freq, 1), 'resting': pct(resting(Pn), 1), 'eig': r(ev, 3), 'memoryYears': r(-1 / np.log(ev[1]), 1), 'reverseDiff': r(np.abs(reverse(Pn) - Pn).max(), 3), 'lags': lags, 'observed': pct([same_record(states, k, wts) for k in lags], 1), 'chain': pct([same_chain(Pn, freq, k) for k in lags], 1), 'chance': pct((freq ** 2).sum(), 1), 'sameNextYear': pct(same_record(states, 1, wts), 1), 'horizonSteps': horizon(Pn)} # (1) the chain of leading doshas, pooled over every place P1, lead_d = describe(leadY, Wn, 4) # (2) the condition of a place: each dosha's strength there with its forty-year mean and straight-line drift taken out, # counted in units of its own usual swing; the state is the dosha that stands highest t = (yrs - yrs.mean())[:, None, None, None] mean = SY.mean(0) slope = (t * (SY - mean)).sum(0) / (t ** 2).sum() resid = SY - mean - slope * t sd = resid.std(0) Wk = Wn * (sd > 0).all(-1) cond = (resid / np.where(sd > 0, sd, 1)).argmax(-1) P2, cond_d = describe(cond, Wk, 3) # the same without taking the drift out, and by the plain ratio to the place's own forty-year measure cond_ratio = (SY / np.where(mean > 0, mean, 1)).argmax(-1) _, ratio_d = describe(cond_ratio, Wk, 3) cond_d['agreesWithRatio'] = pct((Wk[None] * (cond == cond_ratio)).sum() / (Wk.sum() * NY), 0) cond_d['area'] = pct(Wk.sum(), 0) # read only where all three doshas vary from year to year: everywhere but the sea too deep for the column cond_d['areaLand'] = pct((Wk * (near > 0.001)).sum() / (W * (near > 0.001)).sum(), 0) cond_d['areaSea'] = pct((Wk * sea).sum() / (W * sea).sum(), 0) # the same record read apart for land with coast and for the open sea, each against its own chance def kind_record(mask): wk = Wk * mask C1 = chain(cond, wk) f = C1.sum(1) / C1.sum() return {'observed': pct([same_record(cond, k, wk) for k in lags], 1), 'chance': pct((f ** 2).sum(), 1)} cond_d['byKind'] = {'land': kind_record(near > 0.001), 'sea': kind_record(sea)} M['dial'] = {'lead': lead_d, 'condition': cond_d, 'conditionWithDrift': {'lags': lags, 'observed': ratio_d['observed'], 'chance': ratio_d['chance']}, 'back': {'steps': KS, 'lead': {'record': [[r(x, 3) for x in back_record(leadY, k, Wn)] for k in KS], 'chain': [[r(x, 3) for x in np.linalg.matrix_power(reverse(P1), k)] for k in KS]}, 'condition': {'record': [[r(x, 3) for x in back_record(cond, k, Wk)] for k in KS], 'chain': [[r(x, 3) for x in np.linalg.matrix_power(reverse(P2), k)] for k in KS]}}} # (3) month by month: each box against its own calendar month, 60 S to 60 N, where all three doshas vary jj = (lat >= -60) & (lat <= 60) waM, eaM = ground(rainM, P['dry']) Smon = strengths(join(L3, den(sun, 1), waM, eaM)).reshape(NY, 12, NLAT, NLON, 3) clim = Smon.mean(0) ty = (yrs - yrs.mean())[:, None, None, None, None] drift = (ty * (Smon - clim[None])).sum(0) / (ty ** 2).sum() anm = (Smon - clim[None] - drift[None] * ty)[:, :, jj] sdm = anm.std(0) # 12, lat, lon, 3 okm = np.tile((sdm > 1e-9).all(-1)[None], (NY, 1, 1, 1)).reshape(NY * 12, int(jj.sum()), NLON) condm = (anm / np.where(sdm > 1e-9, sdm, 1)[None]).argmax(-1).reshape(NY * 12, int(jj.sum()), NLON) Wj = W[jj] mlags = [1, 2, 3, 4, 6, 9, 12, 18, 24] C3 = chain(condm, Wj, 1, okm) P3, f3 = rows(C3), C3.sum(1) / C3.sum() M['dial']['month'] = {'P': [r(x, 3) for x in P3], 'frequencies': pct(f3, 1), 'lags': mlags, 'observed': pct([same_record(condm, k, Wj, okm) for k in mlags], 1), 'chain': pct([same_chain(P3, f3, k) for k in mlags], 1), 'chance': pct((f3 ** 2).sum(), 1), 'read': pct((okm * Wj[None]).sum() / (NY * 12 * Wj.sum()), 0)} # ---------- for the page ---------- obs = {'lat': [int(x) for x in lat], 'lon': [int(x) for x in lon], 'weights': [round(float(x), 6) for x in w / w.sum() / NLON], 'kind': ''.join('L' if land[j, i] else 'S' if sea[j, i] else 'C' for j in range(NLAT) for i in range(NLON)), 'names': NAMES, 'years': [Y0, Y1], 'now': list(NOW), 'scale': 10000, 'column': COLUMN, 'nowShares': b64u16(E), 'planet': [round(float(x), 5) for x in gE], 'alt': dict({nm: {'nowShares': b64u16(shares(*NOW, G=G, column=col)), 'planet': [round(float(x), 5) for x in gshares(*NOW, G=G, column=col)]} for nm, G, col in (('half', 0.5, None), ('double', 2.0, None), ('surface', 1.0, 0.0), ('column1000', 1.0, 1000.0), ('column6000', 1.0, 6000.0))}, warm={'nowShares': b64u16(Ewarm), 'planet': [round(float(x), 5) for x in gEwarm]}, heat={'nowShares': b64u16(Eheat), 'planet': [round(float(x), 5) for x in gEheat]}, move={'nowShares': b64u16(Emove), 'planet': [round(float(x), 5) for x in gEmove]}), 'landShare': b64u8(near, 250), 'wet': b64u8(wetN, 250), 'monthScale': 250, 'monthShares': b64u8(Em12), 'yearScale': 250, 'yearShares': b64u8(EY), 'globalYears': [[round(float(v), 4) for v in x] for x in gY]} os.makedirs(out_dir, exist_ok=True) json.dump(obs, open(os.path.join(out_dir, 'obs.json'), 'w'), separators=(',', ':')) def clean(x): # a correlation or ratio that cannot be formed (nothing to divide by) is written as null if isinstance(x, dict): return {k: clean(v) for k, v in x.items()} if isinstance(x, list): return [clean(v) for v in x] return None if isinstance(x, float) and x != x else x json.dump(clean(M), open(os.path.join(out_dir, 'obs_measured.json'), 'w'), indent=1) print('wrote', os.path.join(out_dir, 'obs.json'), os.path.getsize(os.path.join(out_dir, 'obs.json')), 'bytes') if __name__ == '__main__': main()