def add_census_data(output):
    # Estimate ANC/SMD Census population data by averaging across the tracts that
    # each ANC and SMD intersects with. Tracts appear to be the lowest level that
    # interesting data is uniformly available from the 2011 American Community Survey.
    # We can get more precise information on population from block groups though.

    # Pre-load the relevant data points from the Census API.

    census_fields = {
            "2010/sf1": [
                    ("P0010001", int, "count"), # Total Population
                    ("P0180001", int, "count"), # Households
                    ("P0180002", int, "count"), # Family Households
                    ("H0040001", int, "count"), # Occupied Housing Units
                    ("H0050001", int, "count"), # Vacant Housing Units
            ],

            # http://www.census.gov/developers/data/acs_5yr_2011_var.xml
            "2011/acs5": [
                    ("B01003_001E", int, "count"), # Total Population
                    ("B01002_001E", float, "median"), # Median Age
                    ("B07001_001E", int, "count"), # Geographic Mobility in the Past Year Total Population
                    ("B07001_017E", int, "count"), # .... Same House One Year Ago
                    ("B07001_033E", int, "count"), # .... Moved Within Same County
                    ("B07001_065E", int, "count"), # .... Moved From Different State
                    ("B07001_081E", int, "count"), # .... Moved From Abroad
                    ("B19019_001E", float, "median"), # Median Household Income
                    ("B16002_001E", int, "count"), # No One Age 14+ in Household Speaks English
            ],
    }

    census_data = { }
    for division in ("tract", "blockgroup"):
        dat = json.load(open("data/smd-%s.json" % division))
        for intsect in dat:
            #if intsect["smd"]["id"] != "1A02": continue

            # get the ID of this tract or blockgroup
            id = intsect[division]["id"]
            if id in census_data: continue
            census_data[id] = { }

            # query Census API
            state_fips, county_fips, tract, bg_num = id[0:2], id[2:5], id[5:11], (id+"?")[11:12]
            for census_table in census_fields:
                field_list = ",".join(f[0] for f in census_fields[census_table])
                if division == "tract":
                    url = "http://api.census.gov/data/%s?key=%s&get=%s&in=state:%s+county:%s&for=tract:%s" \
                        % (census_table, census_api_key, field_list, state_fips, county_fips, tract)
                else:
                    url = "http://api.census.gov/data/%s?key=%s&get=%s&in=state:%s+county:%s+tract:%s&for=blockgroup:%s" \
                        % (census_table, census_api_key, field_list, state_fips, county_fips, tract, bg_num)
                data_from_api = json.load(urlopen(url))
                for i, (fieldname, datatype, sum_mode) in enumerate(census_fields[census_table]):
                    v = data_from_api[1][i]
                    if v is not None: v = datatype(v)
                    census_data[id][fieldname] = v

    # Clear out existing data.
    for ward in output.values():
        for anc in ward["ancs"].values():
            anc["census"] = { "by-blockgroup": { }, "by-tract": { } }
            for smd in anc["smds"].values():
                smd["census"] = { "by-blockgroup": { }, "by-tract": { } }

    # Estimate values for the ANCs and SMDs as a whole.
    for division1 in ("smd", "anc"):
        for division2 in ("tract", "blockgroup"):
            dat = json.load(open("data/%s-%s.json" % (division1, division2)))
            for ix in dat:
                # Here's an intersection between an ANC/SMD and a blockgroup/tract.

                # Get the dict that represents the ANC or SMD's census info.
                ward = ix[division1]["id"][0]
                anc = ix[division1]["id"][1]
                if division1 == "anc":
                    feature = output[ward]["ancs"][anc]
                else:
                    smd = ix[division1]["id"][2:]
                    feature = output[ward]["ancs"][anc]["smds"][smd]
                feature = feature["census"]["by-" + division2]

                # Get the id of the blockgroup or tract.
                census_id = ix[division2]["id"]
                #if census_id not in census_data: continue # only for testing

                # Estimate each field for the ANC/SND as a whole.
                for census_table in census_fields:
                    for fieldname, datatype, summode in census_fields[census_table]:
                        if summode == "count":
                            # This value may not available from the Census (some blockgroup ACS data) so
                            # mark the feature as not-computable.
                            if census_data[census_id][fieldname] is None:
                                feature[fieldname] = "missing-data"
                                continue
                            if feature.get(fieldname) == "missing-data":
                                # Already marked this field is non-computable because data is missing?
                                continue

                            # This field is a count, so we estimate the whole by adding the
                            # field values weighted by the ratio of the intersection area to
                            # the area of the blockgroup or tract.
                            feature.setdefault(fieldname, { "value": 0, "type": summode })
                            feature[fieldname]["value"] += census_data[census_id][fieldname] * ix[division2]["ratio"]
                        elif summode == "median":
                            # For medians, we take the average across all of the intersections
                            # weighted by the population in each intersction.
                            feature.setdefault(fieldname, { "value": 0.0, "type": summode, "weight": 0.0, "missing_weight": 0.0 })
                            w = census_data[census_id]["B01003_001E"] * ix[division2]["ratio"]

                            # This value may not available from the Census (some blockgroup ACS data). It seems
                            # like even tracts are missing some of these values, but we want to estimate as
                            # best as possible so we'll estimate from those tracts where values exist, and we'll
                            # record the proportion of the population in this ANC/SMD that did not contribute
                            # to the estimate.
                            if census_data[census_id][fieldname] is None:
                                feature[fieldname]["missing_weight"] += w
                                continue

                            feature[fieldname]["value"] += census_data[census_id][fieldname] * w
                            feature[fieldname]["weight"] += w

    # Finish up.
    def clean_up(f):
        # Prefer blockgroup statistics if available.
        # Round some values.
        # For medians, divide by the total weight.
        for division in ("tract", "blockgroup"):
            for k, v in list(f["by-" + division].items()):
                if v == "missing-data": continue

                if v["type"] == "count":
                    f[k] = v
                    f[k]["source"] = division
                    v["value"] = int(round(v["value"]))
                elif v["type"] == "median" and v["weight"] > 0:
                    # Skip this if we're missing more data as a blockgroup than as a tract,
                    # or if we're just missing a lot of data (more than 33%).
                    v["missing_pct"] = v["missing_weight"] / (v["weight"] + v["missing_weight"])
                    if v["missing_pct"] > .33 and (k in f and f[k]["missing_pct"] < v["missing_pct"]): continue

                    f[k] = v
                    f[k]["source"] = division
                    v["value"] = round(v["value"] / v["weight"])
                    del v["weight"]
                    del v["missing_weight"]
                del v["type"]
            del f["by-" + division]

    for ward in output.values():
        for anc in ward["ancs"].values():
            clean_up(anc["census"])
            for smd in anc["smds"].values():
                clean_up(smd["census"])
dataset/ETHPy150Open codefordc/ancfinder/scripts/update_anc_database.py/add_census_data/operator `or` -> `and`