Os seis diamantes, lapidados por inteiro

Tudo começa numa escolha: sobre QUE Terra se mede. O app usa o elipsoide GRS80 — semieixo maior a = 6.378.137 m e achatamento f = 1/298,257222101, do qual sai a primeira excentricidade ao quadrado e² = 2f − f². Cada vértice chega em coordenadas geográficas (latitude φ, longitude λ, altura h) no referencial SIRGAS2000.

O primeiro salto é geográfico → geocêntrico (ECEF, o sistema cartesiano com origem no centro de massa da Terra). Usando a grande normal N(φ) = a/√(1 − e²·sin²φ), tem-se X = (N + h)·cosφ·cosλ, Y = (N + h)·cosφ·sinλ e Z = (N·(1 − e²) + h)·sinφ. É a passagem exata de um ponto sobre o elipsoide para (X, Y, Z) em metros.

O segundo salto é a rotação para o plano topocêntrico local ENU (Leste-Norte-Cima), tangente à Terra na origem do polígono — o primeiro vértice (φ₀, λ₀). Para cada vértice calcula-se o vetor ΔR = R − R₀ e projeta-se: e = −sinλ₀·ΔX + cosλ₀·ΔY; n = −sinφ₀·cosλ₀·ΔX − sinφ₀·sinλ₀·ΔY + cosφ₀·ΔZ. Sobra um mapa plano fiel na vizinhança do lote, em que Leste e Norte são metros de verdade.

Com os vértices já no plano (E, N), a área fecha pela fórmula de Gauss (shoelace): A = ½·|Σᵢ (Eᵢ·Nᵢ₊₁ − Eᵢ₊₁·Nᵢ)|, e o perímetro é a soma das distâncias euclidianas entre vértices consecutivos. Quando um lado é curvo (bulge, na convenção DXF), a mesma conta ganha o termo do segmento circular, para a área e o perímetro respeitarem o arco em vez da corda.

A pergunta óbvia é: por que não pegar a área direto do UTM, que já vem em metros? Porque o UTM é uma projeção conforme (preserva ângulos, não áreas) e distorce por um fator de escala k que cresce à medida que se afasta do meridiano central do fuso. A área na carta escala com k², então longe do meridiano ela "incha" — pode passar de meio por cento, o suficiente para uma gleba fechar hectares a mais do que tem no chão.

O SIGEF não aceita essa área distorcida. A norma técnica do INCRA para georreferenciamento de imóveis rurais exige a área calculada no plano topocêntrico local, com precisão posicional certificada (os vértices têm de respeitar tolerâncias em centímetros conforme o tipo de limite). É essa área, e não a do UTM, que vai para o memorial descritivo, do memorial para a matrícula, e da matrícula para a vida jurídica do imóvel — confrontações, retificações, desmembramentos.

Exemplo fictício para fixar: uma "Gleba Modelo" de 12,3456 ha medida no SGL pode aparecer como 12,41 ha se alguém teimar em somar a área direto do UTM a 250 km do meridiano central. Setenta ares de diferença que não existem no terreno — nascidos só da projeção errada. O SGL é o que mantém a peça honesta.

A intuição do plano local é a de quem estende um lençol sobre um morro: bem no ponto onde ele encosta, o lençol e o chão coincidem; só longe dali é que o lençol descola. Medir a área "no lençol" tangente à origem é medir quase sem distorção, porque um lote — mesmo grande — é minúsculo perto do raio da Terra.

A fórmula da área leva o nome de Carl Friedrich Gauss e é conhecida como "fórmula do agrimensor" ou shoelace ("cadarço"), pelo padrão cruzado dos produtos Eᵢ·Nᵢ₊₁. O elipsoide GRS80 foi definido pela União Internacional de Geodésia e Geofísica (IUGG) por volta de 1980 e virou a referência geométrica moderna; o SIRGAS2000 é a materialização sul-americana amarrada ao referencial terrestre internacional (ITRF), que o Brasil oficializou pelo IBGE.

Em sgl.ts, geodeticToECEF faz o primeiro salto (geográfico → geocêntrico) e projetarPlanoLocal faz o segundo (geocêntrico → plano ENU com origem no primeiro vértice); a área fecha por Gauss, com suporte a divisas curvas. O módulo é a fundação de toda peça: é a mesma medida que abastece o memorial, a planilha e a planta.

O código é defensivo de propósito: uma latitude/longitude inválida (NaN) num vértice projeta NaN naquele ponto em vez de cair num fallback de 0°, que jogaria o ponto a milhares de quilômetros e produziria uma área finita e plausível a partir de dado quebrado. A régua da casa é dura: área impossível de calcular tem de aparecer como impossível, nunca virar número certificável calado.

A triangulação de Delaunay é definida por uma propriedade elegante: o circuncírculo de cada triângulo não contém nenhum outro ponto do conjunto — a "propriedade do círculo vazio". Uma consequência famosa é que, entre todas as triangulações possíveis dos mesmos pontos, a de Delaunay MAXIMIZA o menor ângulo. Ou seja, ela evita os triângulos finos e "aguilhoados" e prefere os mais gordos e regulares, o que faz a superfície interpolada ficar bem-comportada.

O teste que decide se um ponto está dentro do circuncírculo (o in-circle test) é o sinal de um determinante 4×4 nas coordenadas dos quatro pontos — uma conta puramente algébrica, exata. Trocar arestas até nenhum teste falhar (edge flipping) é uma das formas clássicas de construir a malha; os melhores algoritmos chegam a O(n·log n).

O diagrama de Voronói é o dual geométrico da Delaunay. A célula de um ponto pᵢ é a região do plano mais próxima dele do que de qualquer outro: { x : ‖x − pᵢ‖ ≤ ‖x − pⱼ‖ para todo j }. Cada fronteira entre duas células é um trecho da mediatriz (o lugar dos pontos equidistantes) do par — e ligar os pares vizinhos de Voronói dá exatamente as arestas de Delaunay.

Um levantamento entrega pontos soltos: vértices, marcos, pontos de apoio, cotas espalhadas. Para desenhar relevo (curvas de nível), interpolar altitude ou raciocinar sobre "de qual marco este trecho está mais perto", é preciso transformar essa nuvem em estrutura. A Delaunay dá a malha triangular (a base de um TIN — rede irregular de triângulos), e o Voronói dá as regiões de influência/proximidade.

A qualidade dos triângulos não é estética: triângulos finos amplificam erro numérico na interpolação de cota e produzem curvas de nível "quebradas". Por isso a garantia de Delaunay — maximizar o menor ângulo — é diretamente uma garantia de altimetria mais estável, que é o que sustenta perfis, seções e a leitura do relevo na planta.

A intuição do Voronói é a do "posto de saúde mais próximo": pinte cada ponto do mapa com a cor do posto mais perto, e as manchas de cor são as células. As fronteiras caem exatamente no meio do caminho entre dois postos — as mediatrizes. A Delaunay é o desenho das amizades: ligue dois postos sempre que suas manchas se encostam.

Georgy Voronoy trabalhou o tema no contexto das formas quadráticas, e a construção leva seu nome desde 1908; a ideia aparece antes em Dirichlet (por isso "tesselação de Dirichlet") e até em Descartes. Boris Delaunay, seu aluno, formalizou a triangulação pela esfera vazia. Como curiosidade, Delaunay foi um alpinista célebre — há um pico com o nome dele.

Em voronoi.ts, o app constrói a célula de cada gerador por interseção de semiplanos recortada numa moldura retangular: a célula do ponto i é a moldura cortada pela mediatriz contra cada outro ponto, ficando sempre o lado mais próximo de i (recorte de Sutherland–Hodgman, exato por construção). Serve para "qual ponto de apoio está mais perto", divisão por proximidade e área de influência por marco.

O módulo é puro e desconfiado: entrada inválida devolve um erro declarado, e ponto duplicado é recusado com o índice dos dois envolvidos — porque a célula de um ponto duplicado é vazia e sumiria calada. Quando a Delaunay entrar de vez na árvore, ela acelera esse dual sem mudar a API.

A krigagem trata a altitude como uma variável regionalizada: um campo em que pontos próximos se parecem e pontos distantes se soltam. Essa "parecença" é medida pelo semivariograma γ(h) = ½·E[(Z(x+h) − Z(x))²] — metade da variância esperada da diferença entre dois pontos separados por uma distância h. Ele sobe da pepita (nugget, o ruído em h→0) até o patamar (sill), e a distância em que estabiliza é o alcance (range). O app modela γ(h) nas três formas clássicas: esférica, exponencial e gaussiana.

O estimador é uma média ponderada dos vizinhos, Ẑ(x₀) = Σ wᵢ·Z(xᵢ). O que faz dele o BLUP (Best Linear Unbiased Predictor) é impor duas coisas: não-viés, forçando Σ wᵢ = 1 (por isso "ordinária": a média é desconhecida), e mínima variância do erro. Isso vira um sistema linear aumentado com um multiplicador de Lagrange μ: para cada vizinho i, Σⱼ wⱼ·γ(dᵢⱼ) + μ = γ(dᵢ₀); e a linha extra Σ wⱼ = 1.

O bônus — e o que separa a krigagem de uma média ponderada qualquer — é que ela devolve a variância do erro de estimativa: σ²ₒₖ = Σ wᵢ·γ(dᵢ₀) + μ. É o número que transforma "a cota aqui é 812,40 m" em "é 812,40 m, e eis o desvio-padrão dessa afirmação". Onde há amostra vizinha, a incerteza cai; num vazio, ela cresce e avisa.

Altimetria vive de interpolar: raramente se mede a cota exatamente onde ela é preciso. O IDW (média por inverso da distância) dá um número, mas cala a dúvida. A krigagem dá o número E a barra de erro — e isso muda a decisão de campo: uma cota estimada com desvio de poucos centímetros num terreno bem amostrado é confiável; a mesma estimativa num rincão sem pontos vizinhos deve acender a luz amarela antes de virar curva de nível na planta.

Por ser um interpolador exato, a krigagem respeita as medições: pedir a cota exatamente sobre uma amostra devolve a própria cota, com variância zero. Isso é o que se espera de uma peça técnica — o modelo não pode "discordar" de um ponto que foi de fato levantado no campo.

A intuição: para adivinhar a altura de um lugar, pese os vizinhos — mas não só pela distância, e sim pela distância comparada ao alcance em que o terreno "esquece" de si mesmo. Se dois vizinhos estão ambos além do alcance, sua opinião vale pouco; se estão dentro dele, valem muito. O variograma é o mapa dessa memória do terreno.

A história é uma bonita passagem da mina para a matemática. Nas minas de ouro de Witwatersrand, na África do Sul, Danie Krige precisava estimar o teor do minério onde a sonda não chegava, e achou um jeito prático de ponderar as amostras. Georges Matheron, na França, pegou aquilo, deu-lhe fundamento probabilístico (a Teoria das Variáveis Regionalizadas), criou a geoestatística e chamou o método de "krigagem", em homenagem a Krige — um método de mineração que virou matemática de superfície.

Em kriging.ts, o app faz krigagem ordinária de cota: monta o sistema pelo semivariograma escolhido (esférico, exponencial ou gaussiano), resolve com o multiplicador de Lagrange por eliminação de Gauss com pivoteamento, e devolve o valor estimado, a variância do erro e os pesos por amostra (que somam 1, úteis para depurar quem puxou a estimativa).

A honestidade está nos casos-limite: alvo coincidente com uma amostra devolve a própria cota com variância zero (interpolador exato); sistema singular — amostras duplicadas ou degeneradas — devolve null, em vez de um número chutado em silêncio. É o degrau acima do IDW de hoje: o mesmo palpite, agora com o tamanho da dúvida grudado nele.

Um levantamento mede coisas demais de propósito: distâncias, ângulos e desníveis redundantes. Como toda medida traz erro, essas observações não fecham exatamente — sobra um resíduo. O método dos mínimos quadrados escolhe os parâmetros (as coordenadas ajustadas) que minimizam a soma ponderada dos quadrados dos resíduos: min vᵀ·P·v, onde v = A·x − w são os resíduos, A é a matriz de desenho (como cada observação depende das incógnitas), w é o vetor de fechamento e P é a matriz de pesos — o inverso da covariância, isto é, medida mais precisa pesa mais.

Derivar e igualar a zero leva às equações normais: (AᵀPA)·x̂ = AᵀPw, cuja solução é x̂ = (AᵀPA)⁻¹·AᵀPw. Os resíduos ajustados v saem no mesmo passo, e a qualidade global vem do desvio-padrão da unidade de peso a posteriori, σ₀² = (vᵀPv)/(n − u), onde n − u são os graus de liberdade (observações menos incógnitas). Sem redundância (n = u), a rede ajusta mas não sobra informação para julgar a qualidade — e σ₀ fica indefinido, honestamente.

Quando as observações são não lineares nas coordenadas (ângulos, redes 2D completas), lineariza-se em torno de valores aproximados e itera-se: é o Gauss–Newton, resolvendo uma sequência de sistemas lineares como o acima até os resíduos pararem de encolher.

Sem ajuste, o erro de fechamento de uma poligonal é empurrado para o último lado, ou espalhado por regra de bolso (a "regra do compasso"). Funciona, mas é cego: não usa a precisão de cada medida e não diz onde está o problema. O ajuste por mínimos quadrados distribui o erro pela rede inteira na proporção do peso de cada observação — e, de brinde, entrega os resíduos, que denunciam qual medida destoou. É o "santo graal" do trabalho de campo: em vez de esconder o erro, mapeá-lo.

Isso é precisão certificável: os resíduos e o σ₀ dão a base estatística para afirmar que os vértices respeitam as tolerâncias exigidas. Uma medida com resíduo grande vira suspeita objetiva de erro grosseiro (um alvo mal apontado, uma anotação trocada), a ser refeita — antes de a peça ir para o INCRA, não depois.

A intuição: imagine as medidas puxando as coordenadas por elásticos, cada elástico mais duro quanto mais precisa a medida. O ponto de equilíbrio — onde a energia total dos elásticos é mínima — é a solução de mínimos quadrados. Elevar ao quadrado (em vez de somar os módulos) é o que torna esse equilíbrio único, suave e resolvível por álgebra linear.

A autoria é uma das disputas mais famosas da matemática. Legendre PUBLICOU o método primeiro, em 1805, com o nome que usamos até hoje. Gauss respondeu que já o usava desde 1795 e o ligou à distribuição normal e à sua célebre previsão de onde reapareceria Ceres em 1801. Os dois têm razão em partes diferentes; a honestidade histórica é dizer exatamente isso — Legendre publicou antes, Gauss aprofundou o fundamento probabilístico.

Em ajusteRede.ts vive o núcleo linear ponderado, genérico e puro: recebe observações e injunções e devolve os parâmetros ajustados, os resíduos por observação e o σ₀, resolvendo as normais por eliminação de Gauss com pivoteamento parcial. Sobre esse núcleo estão as duas aplicações do caso simples: o NIVELAMENTO (a rede 1D clássica de desníveis) e o FECHO de poligonal (dE/dN por lado).

A filosofia da casa aparece de novo: funções puras, sem estado nem UI, e entrada inválida devolvendo um erro declarado — nunca um número plausível calado. As redes 2D completas, não lineares, entram por cima deste núcleo com iteração de Gauss–Newton.

O problema inverso pede: dados dois pontos (φ₁, λ₁) e (φ₂, λ₂), qual a menor distância entre eles SOBRE O ELIPSOIDE (a geodésica) e quais os azimutes de saída e de chegada? Numa esfera a resposta seria a trigonometria esférica; no elipsoide, achatado, a geodésica não é um círculo e a conta fica mais fina.

A ideia clássica (Bessel, Helmert) é levar o problema do elipsoide para uma esfera auxiliar, trocando a latitude φ pela latitude reduzida β, com tanβ = (1 − f)·tanφ. Nessa esfera, o comprimento vira uma integral: s/b = ∫ √(1 + k²·sin²σ) dσ, com k² = e′²·cos²α₀, onde α₀ é o azimute em que a geodésica cruzaria o equador e e′² é a segunda excentricidade. Como k² é minúsculo (da ordem de 0,0067), o integrando é suave e a integral se resolve com pouquíssimos nós.

O caso geral se fecha por iteração: um passo de Newton no azimute de partida α₁ até casar a diferença de longitude pedida. Os casos especiais são exatos: dois pontos no equador têm a própria linha do equador como geodésica; mesma longitude cai numa integral ao longo do meridiano.

O memorial descritivo não pede a distância "em linha reta na tela" nem sobre uma esfera aproximada: pede a distância e o azimute geodésicos, sobre o elipsoide de referência. Aproximar a Terra por esfera introduz erro de dezenas de metros em lados longos — inaceitável numa peça em que o vértice tem tolerância de centímetros. É por isso que o azimute do memorial só "fecha" com a conta feita no elipsoide.

E há a robustez: o método antigo de Vincenty pode não convergir quando os dois pontos são quase antipodais (quase do outro lado do mundo). Numa peça oficial, um número errado saindo calado é o pior cenário. Por isso Karney importa duplamente — precisão e a garantia de que, quando não dá para confiar, o método avisa em vez de mentir.

A intuição da geodésica: é o caminho que um fio bem esticado faria sobre a casca da Terra — reto no sentido de "sem virar para lado nenhum", ainda que curvo visto de fora. Sobre uma laranja levemente achatada, esse fio não fecha um círculo perfeito, e a distância depende sutilmente de por onde ele passa; capturar isso é o trabalho das séries e integrais.

Charles Karney é o único vivo desta lista. Em 2013 publicou algoritmos (reunidos na biblioteca GeographicLib) que medem distância e rumo no elipsoide com precisão de nanômetros e — diferente de Vincenty — sem falhar nos casos quase antipodais. É a modernização de uma linhagem que passa por Bessel e Helmert no século XIX.

Em geodesicaKarney.ts, o app resolve o problema inverso na esfera auxiliar de Karney, mas — em vez de reconstruir de memória as séries truncadas do GeographicLib, o que arriscaria um número errado numa peça — resolve os integrais EXATOS por quadratura de Gauss–Legendre de 16 nós. O integrando é tão suave (k² ≤ e′² ≈ 0,0067) que isso já dá muito abaixo do milímetro.

O módulo trata os casos com cuidado: equatorial e meridiano por fórmula fechada, geral por Newton no azimute de partida, e — o ponto sagrado — quase-antípodas que não convergem retornam confiavel:false, para nunca deixar um número não confiável entrar numa peça. Aqui o app usa as constantes do elipsoide WGS84 (a = 6.378.137 m, f = 1/298,257223563), que difere do GRS80 do SGL apenas na 12ª casa do achatamento — irrelevante para distância e azimute.

Há duas Terras de referência ao mesmo tempo. O ELIPSOIDE (GRS80) é uma superfície geométrica lisa — é sobre ele que o GNSS/GPS entrega a altura, chamada altura elipsoidal h. O GEOIDE é uma superfície física: a superfície equipotencial do campo de gravidade que melhor coincide com o nível médio dos mares — é a ela que se refere a altitude ortométrica H, a "altura acima do nível do mar", a que decide para onde a água corre. A separação entre as duas, ponto a ponto, é a ondulação geoidal N.

A relação fundamental é simples e implacável: h = H + N, logo H = h − N. O GPS mede h; o memorial e a hidrologia querem H; então é preciso conhecer N em cada vértice. No Brasil N é sempre negativo, variando tipicamente de cerca de −4 m no Sul a algo perto de −35 m na Amazônia — ou seja, ignorar N erra a altitude em dezenas de metros.

De onde vem N? Da física do campo de gravidade. A fórmula de Bruns, N = T/γ, liga a ondulação ao potencial perturbador T (o quanto a gravidade real difere da do elipsoide) dividido pela gravidade normal γ; a integral de Stokes calcula N a partir das anomalias de gravidade medidas no mundo todo. Na prática, isso já vem embalado em modelos: o MAPGEO do IBGE e o EGM2008 (uma expansão em harmônicos esféricos do campo de gravidade até grau altíssimo).

A altitude entra na área. No SGL, cada vértice é levado ao geocêntrico com sua altura h, então um erro sistemático na altura contamina, ainda que pouco, a geometria projetada. Mais direto: a altimetria da peça — cotas, curvas de nível, perfis — precisa estar em H (ortométrica), porque é ela que tem sentido físico. Uma curva de nível desenhada com h elipsoidal estaria dezenas de metros fora do que a régua de campo e o nível medem.

Um DEM (modelo digital de elevação) como o do Copernicus/SRTM entrega H; o receptor entrega h. Para casar os dois mundos é preciso somar ou subtrair N no lugar certo — e "no lugar certo" muda de município para município. Errar N é um deslocamento silencioso que entra em todas as cotas do imóvel de uma vez.

A intuição: o elipsoide é a Terra idealizada, uma bola lisa de sabonete; o geoide é a Terra que a água conhece, cheia de altos e baixos onde a gravidade é mais forte ou mais fraca (uma montanha densa "puxa" o nível para cima; uma bacia, para baixo). A ondulação N é a distância entre esses dois retratos — e ela é sempre suave, porque a gravidade varia devagar no espaço.

A história atravessa o século XIX. Gauss já falava da "figura matemática da Terra" como a superfície de nível dos oceanos idealizada; foi Johann Listing quem, em 1873, deu-lhe o nome "geoide". George Stokes forneceu, em 1849, a integral que reconstrói o geoide a partir da gravidade, e Bruns, a relação com o potencial. É a geodésia física clássica que hoje mora dentro de um arquivo do app.

Em geoid.ts, o app interpola N para qualquer coordenada do Brasil por IDW (média por inverso da distância, com potência de decaimento) sobre um conjunto de âncoras MAPGEO/EGM distribuídas pelas cinco regiões, e converte nos dois sentidos: elipsoidalParaOrtometrica (H = h − N) e ortometricaParaElipsoidal (h = H + N). Esse N é somado a cada altitude baixada do DEM antes de a cota entrar na área do SGL.

A honestidade aparece no valor de retorno: coordenada inválida devolve NaN, não um antigo "valor médio Brasil" de −15 m que tinha a mesma cara de uma ondulação legítima e entraria calado em todas as cotas. Com NaN, a cota deixa de ser finita, o portão de exportação reconhece a altitude ausente e trava a peça — em vez de deixar sair uma área plausível e errada. Número não calculado não vira número.