Maneira mais precisa de calcular a área de rasters

No meu trabalho diário, sou constantemente solicitado a calcular áreas de conjuntos de dados raster globais em projeção geográfica com resolução de 30 segundos de arco. Esses conjuntos de dados normalmente são o resultado de uma operação Combinar (um exemplo típico é uma classe de vegetação combinada com a camada do país). Para fazer isso, nossa unidade criou um conjunto de dados raster com a área de cada pixel na projeção geográfica a 30 arco-segundos. Com essa grade de área, é executado um estado zonal para somar as áreas de cada classe. Como não tenho certeza de como essa grade de área foi criada, sempre me perguntei se essa abordagem é mais precisa de reprojetar o raster em uma projeção de área igual (a partir de testes simples, os resultados dos dois métodos são semelhantes). Alguém experimentou uma situação semelhante?

Existe uma fórmula exata relativamente simples para a área de qualquer quadrângulo esférico delimitado por paralelos (linhas de latitude) e meridianos (linhas de longitude). Ele pode ser derivado diretamente usando propriedades básicas da elipse (do eixo maior a e do eixo menor b ) que é girada em torno de seu eixo menor para produzir o elipsóide. (A derivação é um bom exercício de cálculo integral, mas acredito que seria de pouco interesse neste site.)

A fórmula é simplificada dividindo o cálculo em etapas básicas.

Primeiro, a distância entre os limites leste e oeste - os meridianos l0 e l1 - é uma fração de um círculo inteiro igual a q = (l1 - l0) / 360 (quando os meridianos são medidos em graus) ou 1 = ( l1 - l0) / (2 * pi) (quando os meridianos são medidos em radianos). Encontre a área de toda a fatia localizada entre os paralelos f0 e f1 e apenas multiplique por q .

Segundo, empregaremos uma fórmula para a área de uma fatia horizontal do elipsóide delimitada pelo Equador (em f0 = 0) e um paralelo na latitude f (= f1). A área da fatia entre duas latitudes f0 ef1 (situada no mesmo hemisfério) será a diferença entre a área maior e a menor.

Finalmente, desde que o modelo seja realmente um elipsóide (e não uma esfera), a área dessa fatia entre o Equador e o paralelo na latitude f é dada por

area(f) = pi * b^2 * (log(zp/zm) / (2*e) + sin(f) / (zp*zm))

onde ae bsão os comprimentos dos eixos maior e menor da elipse geradora, respectivamente,

e = sqrt(1 - (b/a)^2)

é sua excentricidade e

zm = 1 - e*sin(f); zp = 1 + e*sin(f)

(Isso é muito mais simples do que computar com geodésica, que são apenas aproximações aos paralelos de qualquer maneira. Observe o comentário de @cffk sobre uma maneira de calcular de log(zp/zm)maneira a evitar perda de precisão em baixas latitudes.)

area(f) é a área da fatia opaca do equador até a latitude f (cerca de 30 graus norte na ilustração. X e Y são eixos de coordenadas cartesianas geocêntricas mostrados para referência.

Para o elipsóide WGS 84, use os valores constantes

a = 6 378 137 meters,  b =  6 356 752.3142 meters,

implicando

e = 0.08181919084296

(Para um modelo esférico com a = b , a fórmula se torna indefinida. É necessário definir um limite como e -> 0 a partir de cima, o que reduz a fórmula padrão 2 * pi * a^2 * sin(f).)

De acordo com essas fórmulas, um quadrângulo de 30 'por 30' baseado no Equador tem uma área de 3077,2300079129 quilômetros quadrados, enquanto um quadrângulo de 30 'por 30' tocando um poste (que na verdade é apenas um triângulo) tem uma área de apenas 13,6086152 quilômetros.

Como verificação, as fórmulas aplicadas a todas as células de uma grade de 720 por 360 cobrindo a superfície da Terra fornecem uma área total de 4 * pi * (6371.0071809) ^ 2 quilômetros quadrados, indicando que o raio autálico da Terra deve ser 6371.0071809 quilômetros. Isso difere do valor da Wikipedia apenas na última figura significativa (cerca de um décimo de milímetro). (Eu acho que os cálculos da Wikipedia estão um pouco errados :-).

Como verificações adicionais, usei versões dessas fórmulas para reproduzir os Apêndices 4 e 5 em Lev M. Bugayevskiy e John P. Snyder, Map Projections: A Reference Manual (Taylor & Francis, 1995). O Apêndice 4 mostra comprimentos de arco de seções de 30 'de meridianos e paralelos, dados ao medidor mais próximo. Uma verificação pontual dos resultados mostrou concordância perfeita. Depois recriei a tabela com incrementos de 0,0005 ', em vez de 0,5', e integrei numericamente as áreas dos quadrângulos, conforme estimado com esses comprimentos de arco. A área total do elipsóide foi reproduzida com precisão para melhor que oito números significativos. O Apêndice 5 mostra os valores de area(f)para f = 0, 1/2, 1, ..., 90 graus, multiplicados por 1 / (2 * pi). Esses valores são dados ao quilômetro quadrado mais próximo. Uma verificação visual de valores próximos de 0, 45 e 90 graus mostrou concordância perfeita.

Essa fórmula exata pode ser aplicada usando álgebra raster começando com uma grade que fornece as latitudes dos limites superiores de cada célula e outra que indica as latitudes dos limites inferiores. Cada uma delas é, essencialmente, uma grade de coordenadas y. (Em cada caso, você pode querer criar sin(f)e então zme zpcomo resultados intermediários.) Subtraia os dois resultados, pegue o valor absoluto disso e multiplique pela fração q obtida na primeira etapa (igual a 0,5 / 360 = 1/720 para uma largura de célula de 30 pés, por exemplo). Esta será uma grade cujos valores contêm exatamenteáreas de cada célula (até a precisão numérica da própria grade). Apenas certifique-se de expressar as latitudes na forma esperada pela função seno: muitas calculadoras raster fornecerão coordenadas em graus, mas esperam radianos por suas funções trigonométricas!

Para o registro, aqui estão as áreas exatas de células 30 'por 30' no elipsóide WGS 84 do Equador até um polo, em intervalos de 30 ', a 11 figuras (o mesmo número usado para o raio menor b ):

3077.2300079,3077.0019391,3076.5458145,3075.8616605,3074.9495164,3073.8094348,3072.4414813,3070.8457347,3069.0222870,3066.9712434,3064.6927222,3062.1868550,3059.4537865,3056.4936748,3053.3066912,3049.8930202,3046.2528597,3042.3864209,3038.2939285,3033.9756204,3029.4317480,3024.6625762,3019.6683833,3014.4494612,3009.0061153,3003.3386648,2997.4474422,2991.3327939,2984.9950800,2978.4346744,2971.6519646,2964.6473522,2957.4212526,2949.9740951,2942.3063230,2934.4183938,2926.3107788,2917.9839636,2909.4384482,2900.6747464,2891.6933866,2882.4949115,2873.0798782,2863.4488581,2853.6024374,2843.5412166,2833.2658109,2822.7768503,2812.0749792,2801.1608571,2790.0351582,2778.6985716,2767.1518013,2755.3955665,2743.4306011,2731.2576543,2718.8774905,2706.2908892,2693.4986451,2680.5015685,2667.3004848,2653.8962347,2640.2896746,2626.4816763,2612.4731271,2598.2649300,2583.8580035,2569.2532818,2554.4517149,2539.4542684,2524.2619238,2508.8756783,2493.2965451,2477.5255533,2461.5637477,2445.4121891,2429.0719545,2412.5441367,2395.8298444,2378.9302026,2361.8463521,2344.5794500,2327.1306692,2309.5011988,2291.6922441,2273.7050264,2255.5407830,2237.2007674,2218.6862492,2199.9985139,2181.1388633,2162.1086151,2142.9091030,2123.5416769,2104.0077025,2084.3085615,2064.4456516,2044.4203864,2024.2341953,2003.8885234,1983.3848318,1962.7245972,1941.9093120,1920.9404843,1899.8196375,1878.5483108,1857.1280585,1835.5604507,1813.8470724,1791.9895239,1769.9894206,1747.8483931,1725.5680867,1703.1501618,1680.5962932,1657.9081707,1635.0874985,1612.1359952,1589.0553936,1565.8474409,1542.5138984,1519.0565410,1495.4771578,1471.7775513,1447.9595378,1424.0249466,1399.9756206,1375.8134157,1351.5402005,1327.1578567,1302.6682785,1278.0733724,1253.3750574,1228.5752643,1203.6759360,1178.6790272,1153.5865040,1128.4003439,1103.1225355,1077.7550785,1052.2999830,1026.7592702,1001.1349711,975.42912705,949.64378940,923.78101904,897.84288636,871.83147097,845.74886152,819.59715539,793.37845851,767.09488512,740.74855748,714.34160569,687.87616739,661.35438752,634.77841811,608.15041795,581.47255240,554.74699308,527.97591765,501.16150951,474.30595754,447.41145586,420.48020351,393.51440422,366.51626611,339.48800143,312.43182627,285.34996030,258.24462644,231.11805066,203.97246162,176.81009042,149.63317034,122.44393648,95.244625564,68.037475592,40.824725575,13.608615243

Os valores estão em quilômetros quadrados.

Se você deseja aproximar essas áreas ou simplesmente entender melhor seu comportamento, a fórmula se reduz a uma série de potências seguindo esse padrão:

area(f) = 2 * pi * b^2 * z * (1 + (4/3)y + (6/5)y^2 + (8/7)y^3 + ...)

Onde

z = sin(f), y = (e*z)^2.

(Uma fórmula equivalente aparece em Bugayevskiy & Snyder, op. Cit. , Equação (2.1).)

Como e ^ 2 é muito pequeno (em torno de 1/150 para todos os modelos elipsoidais da Terra) e z fica entre 0 e 1, y também é pequeno. Assim, os termos y ^ 2, y ^ 3, ... rapidamente diminuem, adicionando mais de duas casas decimais com precisão a cada termo. Se ignorássemos y completamente, a fórmula seria a da área de uma esfera de raio b . Os termos restantes podem ser entendidos como correção da protuberância equatorial da Terra.

Editar

Algumas questões foram levantadas sobre como o cálculo da distância geodésica da área se compara com essas fórmulas exatas. O método da distância geodésica aproxima cada quadrângulo pela geodésica, em vez de paralelos, que conectam seus cantos horizontalmente e aplica a fórmula euclidiana para um trapézio. Para pequenos quadrângulos, como os quadrantes de 30 ', isso é ligeiramente inclinado e tem uma precisão relativa entre 6 e 10 partes por milhão. Aqui está um gráfico do erro para o WGS 84 (ou qualquer elipsóide de terra razoável, para esse assunto):

Portanto, se (1) você tiver acesso fácil aos cálculos de distância geodésica e (2) puder tolerar erros no nível de ppm, considere usar esses cálculos geodésicos e multiplicar seus resultados por 1,00000791 para corrigir o viés. Para mais duas casas decimais de precisão, subtraia pi / 2 * cos (2f) / 10 ^ 6 do fator de correção: o resultado será preciso dentro de 0,04 ppm.

A resposta à pergunta de radouxju depende da forma do pixel quando projetada no elipsóide. Se o sistema de coordenadas da varredura for longitude e latitude, o pixel é um retângulo de linha de rumba e a resposta do whuber pode ser usada ou, de maneira mais geral, você pode usar a fórmula para um polígono cujas bordas são linhas de rumba. Se o sistema de coordenadas for uma projeção conforme em larga escala (UTM, plano de estado etc.), seria mais preciso aproximar as arestas por geodésica e usar a fórmula para um polígono geodésico. Os polígonos geodésicos são provavelmente os melhores para uso geral, pois, diferentemente dos polígonos da linha de rumba, eles são "bem comportados" próximos aos pólos.

Implementações das fórmulas para polígonos geodésicos e de linha de rumba são fornecidas pela minha biblioteca GeographicLib . A área geodésica está disponível em vários idiomas; a área da linha de rhumb é apenas C ++. Há uma versão online (linha geodésica + rhumb) disponível aqui . A precisão desses cálculos é tipicamente melhor que 0,1 metro quadrado.

Você terá que julgar credibilidade / oficial ... As fórmulas geodésicas são derivadas em A área sob a geodésica (Danielsen, 1989, é necessária assinatura) e Algoritmos para geodésica (Karney, 2013, acesso aberto). As fórmulas de linha de rhumb são fornecidas aqui .

Deparei com essa pergunta ao tentar determinar uma fórmula para a área de um pixel WGS84. Embora a resposta do @ whuber contenha essas informações, ainda havia algum trabalho para obter uma fórmula para a área de um pixel de grau quadrado em uma determinada latitude. Eu incluí uma função Python que escrevi abaixo que abstrai isso em uma única chamada. Embora não responda diretamente à pergunta do pôster sobre a área de uma varredura INTEIRA (embora seja possível somar as áreas de todos os pixels), acho que ainda são informações úteis para alguém que possa estar procurando um cálculo semelhante.

def area_of_pixel(pixel_size, center_lat):
    """Calculate m^2 area of a wgs84 square pixel.

    Adapted from: /gis//a/127327/2397

    Parameters:
        pixel_size (float): length of side of pixel in degrees.
        center_lat (float): latitude of the center of the pixel. Note this
            value +/- half the `pixel-size` must not exceed 90/-90 degrees
            latitude or an invalid area will be calculated.

    Returns:
        Area of square pixel of side length `pixel_size` centered at
        `center_lat` in m^2.

    """
    a = 6378137  # meters
    b = 6356752.3142  # meters
    e = math.sqrt(1 - (b/a)**2)
    area_list = []
    for f in [center_lat+pixel_size/2, center_lat-pixel_size/2]:
        zm = 1 - e*math.sin(math.radians(f))
        zp = 1 + e*math.sin(math.radians(f))
        area_list.append(
            math.pi * b**2 * (
                math.log(zp/zm) / (2*e) +
                math.sin(math.radians(f)) / (zp*zm)))
    return pixel_size / 360. * (area_list[0] - area_list[1])