"""Fase 2: baixa e monta as variáveis municipais de exposição definidas em PLANO-fase2.md. Fontes (IBGE/SIDRA, API apisidra): PAM 1612 e 1613, PEVS 289 e 291 (média 2005–2010), Censo 2022 (6803, 9923), PIB municipal 5938 (2010). Altitude do centroide: Open-Meteo (Copernicus DEM 90 m). Tratamento dos símbolos do IBGE (revisto depois da primeira versão; ver auditoria em fase2_auditoria_ausentes.csv): - "-" é zero absoluto e vira 0; - "..." (não disponível), ".." (não se aplica), "X" (sigilo) e municípios ausentes da resposta ficam AUSENTES; - a média de 2005–2010 usa só os anos com dado; sem nenhum ano com dado, o município fica ausente; - lenha = extração + silvicultura, somadas só nos anos em que as duas existem. A versão antiga (ausente = 0) é mantida como sensibilidade em fase2_exposicoes_ausente_zero.csv. Comparação histórica (não pré-especificada): população rural em 2000 e 2010 (tabela 202) e domicílios com poço ou nascente em 2010 (tabela 3217). Saídas: dados/brutos/fase2/*.json (respostas originais), resultados/fase2_exposicoes.csv, resultados/fase2_exposicoes_ausente_zero.csv e resultados/fase2_auditoria_ausentes.csv Uso: python exposicoes.py """ import json import math import time import numpy as np import pandas as pd import requests from analisar import centroides from config import BRUTOS, RESULTADOS PASTA = BRUTOS / "fase2" ANOS = "2005-2010" CONSULTAS = { "fumo_area": f"t/1612/n6/all/v/216/p/{ANOS}/c81/2703", "mate_cultivo_area": f"t/1613/n6/all/v/216/p/{ANOS}/c82/2729", "cafe_area": f"t/1613/n6/all/v/216/p/{ANOS}/c82/2723", "mate_extracao_t": f"t/289/n6/all/v/144/p/{ANOS}/c193/3406", "lenha_extracao_m3": f"t/289/n6/all/v/144/p/{ANOS}/c193/3434", "lenha_silvicultura_m3": f"t/291/n6/all/v/142/p/{ANOS}/c194/3456", "agua_2022": "t/6803/n6/all/v/381/c1821/72129,72146,72147,72148,72154,72155,72156", "situacao_2022": "t/9923/n6/all/v/93/p/2022/c1/6795,2", "pib_2010": "t/5938/n6/all/v/37/p/2010", "situacao_2000_2010": "t/202/n6/all/v/93/p/2000,2010/c2/0/c1/0,2", "agua_2010": "t/3217/n6/all/v/96/c1/0/c61/0,10971,121290,121296/c67/0/c309/0", } ANOS_LISTA = list(range(2005, 2011)) def coluna_ano(df): return next(c for c in df.columns if c.startswith("D") and c.endswith("C") and set(df[c]) <= {str(a) for a in range(1970, 2030)}) def media_anos(df, indice, ausente_zero): """Média 2005–2010 por município. "-" = 0; demais símbolos e linhas ausentes = ausente (ou 0, na sensibilidade).""" d = df.copy() d["ano"] = d[coluna_ano(d)].astype(int) d["x"] = np.where(d["V"].str.strip() == "-", 0.0, d["valor"]) tab = d.pivot_table(index="cod_ibge", columns="ano", values="x", aggfunc="first", dropna=False) tab = tab.reindex(index=indice, columns=ANOS_LISTA) if ausente_zero: tab = tab.fillna(0) return tab def auditoria(df, nome, indice): d = df.copy() d["ano"] = d[coluna_ano(d)].astype(int) tab = d.pivot_table(index="cod_ibge", columns="ano", values="V", aggfunc="first").reindex(index=indice, columns=ANOS_LISTA) v = tab.stack(future_stack=True) return {"variavel": nome, "celulas": v.size, "numericas": int(pd.to_numeric(v, errors="coerce").notna().sum()), "zero_absoluto_traco": int((v == "-").sum()), "nao_disponivel_3pontos": int((v == "...").sum()), "nao_se_aplica_2pontos": int((v == "..").sum()), "sigilo_X": int((v == "X").sum()), "ausente_da_resposta": int(v.isna().sum()), "municipios_sem_nenhum_ano_com_dado": int(tab.apply(lambda r: pd.to_numeric(r.where(r != "-", "0"), errors="coerce").notna().sum() == 0, axis=1).sum())} def baixar(nome, consulta): destino = PASTA / f"{nome}.json" if not destino.exists(): for tentativa in range(4): try: r = requests.get(f"https://apisidra.ibge.gov.br/values/{consulta}", timeout=600) r.raise_for_status() destino.write_text(r.text) break except requests.RequestException: time.sleep(10 * (tentativa + 1)) else: raise RuntimeError(f"falhou: {nome}") linhas = json.loads(destino.read_text())[1:] df = pd.DataFrame(linhas) df["cod_ibge"] = df["D1C"].astype(int) df["valor"] = pd.to_numeric(df["V"], errors="coerce") return df def coluna_categoria(df, codigo): """Nome da coluna D*C que contém o código de categoria (a posição varia entre tabelas).""" return next(c for c in df.columns if c.startswith("D") and c.endswith("C") and codigo in set(df[c])) def area_km2(geom): """Área aproximada do polígono (projeção equirretangular local).""" total = 0.0 polys = geom["coordinates"] if geom["type"] == "MultiPolygon" else [geom["coordinates"]] for poly in polys: for k, anel in enumerate(poly): lat0 = math.radians(sum(p[1] for p in anel) / len(anel)) xs = [math.radians(p[0]) * math.cos(lat0) * 6371 for p in anel] ys = [math.radians(p[1]) * 6371 for p in anel] a = 0.5 * abs(sum(xs[i] * ys[i + 1] - xs[i + 1] * ys[i] for i in range(len(anel) - 1))) total += a if k == 0 else -a return total def altitudes(cent): """Altitude do centroide, em lotes de 100; o progresso fica salvo em altitude.json.""" destino = PASTA / "altitude.json" saida = {int(k): v for k, v in json.loads(destino.read_text()).items()} if destino.exists() else {} faltam = [c for c in sorted(cent) if c not in saida] for i in range(0, len(faltam), 100): lote = faltam[i:i + 100] params = {"latitude": ",".join(f"{cent[c][0]:.4f}" for c in lote), "longitude": ",".join(f"{cent[c][1]:.4f}" for c in lote)} for tentativa in range(6): try: r = requests.get("https://api.open-meteo.com/v1/elevation", params=params, timeout=90) if r.status_code == 200: break except requests.RequestException: pass time.sleep(10 * (tentativa + 1)) else: raise RuntimeError("altitude: API indisponível") saida.update(dict(zip(lote, r.json()["elevation"]))) destino.write_text(json.dumps(saida)) time.sleep(1) return saida def main(): PASTA.mkdir(parents=True, exist_ok=True) geo = json.loads((BRUTOS / "malha_municipios.geojson").read_text()) base = pd.DataFrame({"cod_ibge": [int(f["properties"]["codarea"]) for f in geo["features"]], "area_km2": [area_km2(f["geometry"]) for f in geo["features"]]}).set_index("cod_ibge") p10 = pd.read_csv(BRUTOS / "populacao_censo2010.csv").groupby("cod_ibge")["pop"].sum() p22 = pd.read_csv(BRUTOS / "populacao_censo2022.csv").groupby("cod_ibge")["pop"].sum() base["pop_2010"] = p10.reindex(base.index) base["pop_2022"] = p22.reindex(base.index) pop_ref = base["pop_2010"].fillna(base["pop_2022"]) # municípios criados depois de 2010 usam 2022 hectares = base["area_km2"] * 100 brutos = {nome: baixar(nome, CONSULTAS[nome]) for nome in ["fumo_area", "mate_cultivo_area", "cafe_area", "mate_extracao_t", "lenha_extracao_m3", "lenha_silvicultura_m3"]} pd.DataFrame([auditoria(df, nome, base.index) for nome, df in brutos.items()]).to_csv( RESULTADOS / "fase2_auditoria_ausentes.csv", index=False) versoes = {} for rotulo, zero in (("estrita", False), ("ausente_zero", True)): t = {nome: media_anos(df, base.index, zero) for nome, df in brutos.items()} v = base[[]].copy() v["fumo"] = 100 * t["fumo_area"].mean(axis=1) / hectares v["mate_cultivo"] = 100 * t["mate_cultivo_area"].mean(axis=1) / hectares v["cafe"] = 100 * t["cafe_area"].mean(axis=1) / hectares v["mate_extracao"] = 1000 * t["mate_extracao_t"].mean(axis=1) / pop_ref v["lenha"] = (t["lenha_extracao_m3"] + t["lenha_silvicultura_m3"]).mean(axis=1) / pop_ref v["lenha_extracao"] = t["lenha_extracao_m3"].mean(axis=1) / pop_ref # sensibilidade: só extração versoes[rotulo] = v print(f" {rotulo}: ausentes por variável", v.isna().sum().to_dict()) base = base.join(versoes["estrita"]) agua = baixar("agua_2022", CONSULTAS["agua_2022"]) agua = agua.pivot_table(index="cod_ibge", columns=coluna_categoria(agua, "72129"), values="valor", aggfunc="sum") total = agua["72129"] base["agua_poco"] = (100 * agua[["72146", "72147", "72154", "72155"]].sum(axis=1) / total).reindex(base.index) base["agua_nascente"] = (100 * agua[["72148", "72156"]].sum(axis=1) / total).reindex(base.index) sit = baixar("situacao_2022", CONSULTAS["situacao_2022"]) sit = sit.pivot_table(index="cod_ibge", columns=coluna_categoria(sit, "6795"), values="valor", aggfunc="sum") base["rural"] = (100 * sit["2"] / sit["6795"]).reindex(base.index) pib = baixar("pib_2010", CONSULTAS["pib_2010"]).set_index("cod_ibge")["valor"] base["pib_pc"] = np.log(1000 * pib.reindex(base.index) / pop_ref) sit0 = baixar("situacao_2000_2010", CONSULTAS["situacao_2000_2010"]) sit0["ano"] = sit0[coluna_ano(sit0)] cat = coluna_categoria(sit0, "2") for ano in ("2000", "2010"): x = sit0[sit0["ano"] == ano].pivot_table(index="cod_ibge", columns=cat, values="valor", aggfunc="sum") base[f"rural_{ano}"] = (100 * x["2"] / x["0"]).reindex(base.index) ag10 = baixar("agua_2010", CONSULTAS["agua_2010"]) ag10 = ag10.pivot_table(index="cod_ibge", columns=coluna_categoria(ag10, "10971"), values="valor", aggfunc="sum") base["agua_poco_nascente_2010"] = (100 * ag10[["10971", "121290", "121296"]].sum(axis=1, min_count=1) / ag10["0"]).reindex(base.index) cent = centroides() alt = altitudes(cent) base["altitude"] = pd.Series(alt).reindex(base.index) base.to_csv(RESULTADOS / "fase2_exposicoes.csv", float_format="%.6g") zero = base.copy() for c in versoes["ausente_zero"].columns: zero[c] = versoes["ausente_zero"][c] zero.to_csv(RESULTADOS / "fase2_exposicoes_ausente_zero.csv", float_format="%.6g") print(pd.read_csv(RESULTADOS / "fase2_auditoria_ausentes.csv").to_string(index=False)) pd.set_option("display.width", 200) print(base.describe().T.round(3).to_string()) print("Ausentes:", base.isna().sum()[base.isna().sum() > 0].to_dict()) if __name__ == "__main__": main()