Un moyen plus précis de calculer la surface des rasters

Dans mon travail quotidien, on me demande constamment de calculer des zones de jeux de données raster globales en projection géographique à une résolution de 30 secondes d'arc. Ces ensembles de données sont normalement le résultat d'une opération de combinaison (un exemple typique est une classe de végétation combinée avec une couche de pays). Pour ce faire, notre unité a créé un jeu de données raster avec l'aire de chaque pixel en projection géographique à 30 secondes d'arc. Avec cette grille de surfaces, un zonalstat est effectué pour additionner les surfaces pour chaque classe. Comme je ne sais pas comment cette grille de surface a été créée, je me suis toujours demandé si cette approche était plus précise que de simplement reprojeter le raster dans une projection à surface égale (à partir de tests simples, les résultats des deux méthodes sont similaires). Quelqu'un a-t-il vécu une situation similaire?

Il existe une formule exacte relativement simple pour l'aire de tout quadrilatère sphérique délimité par des parallèles (lignes de latitude) et des méridiens (lignes de longitude). Il peut être dérivé directement en utilisant les propriétés de base de l'ellipse (des axes majeurs a et secondaire b ) qui est tournée autour de son petit axe pour produire l'ellipsoïde. (La dérivation fait un bel exercice de calcul intégral, mais je pense que ce serait peu intéressant sur ce site.)

La formule est simplifiée en décomposant le calcul en étapes de base.

Premièrement, la distance entre les limites est et ouest - les méridiens l0 et l1 - est une fraction d'un cercle entier égale à q = (l1 - l0) / 360 (lorsque les méridiens sont mesurés en degrés) ou 1 = ( l1 - l0) / (2 * pi) (lorsque les méridiens sont mesurés en radians). Trouvez l'aire de la tranche entière située entre les parallèles f0 et f1 et multipliez-la par q .

Deuxièmement, nous utiliserons une formule pour l'aire d'une tranche horizontale de l'ellipsoïde délimitée par l'équateur (à f0 = 0) et d'un parallèle à la latitude f (= f1). La zone de la tranche entre deux latitudes f0 et f1 (se trouvant sur le même hémisphère) sera la différence entre la plus grande et la plus petite zone.

Enfin, à condition que le modèle soit vraiment un ellipsoïde (et non une sphère), l'aire d'une telle tranche entre l'équateur et le parallèle à la latitude f est donnée par

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

où aet bsont les longueurs des axes majeur et mineur de l'ellipse génératrice, respectivement,

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

est son excentricité, et

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

(C'est beaucoup plus simple que de calculer avec des géodésiques, qui ne sont de toute façon que des approximations des parallèles. Veuillez noter le commentaire de @cffk concernant un moyen de calculer log(zp/zm)d'une manière qui évite la perte de précision aux basses latitudes.)

area(f) est l'aire de la tranche opaque de l'équateur jusqu'à la latitude f (environ 30 degrés nord dans l'illustration. X et Y sont des axes de coordonnées cartésiennes géocentriques indiqués à titre de référence.

Pour l'ellipsoïde WGS 84, utilisez les valeurs constantes

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

entraînant

e = 0.08181919084296

(Pour un modèle sphérique avec a = b , la formule devient indéfinie. Vous devez prendre une limite e -> 0 par le haut, qui se réduit ensuite à la formule standard 2 * pi * a^2 * sin(f).)

Selon ces formules, un quadrangle de 30 'par 30' basé sur l'équateur a une superficie de 3077,2300079129 kilomètres carrés, tandis qu'un quadrangle de 30 'par 30' touchant un poteau (qui n'est en fait qu'un triangle) a une superficie de seulement 13,6086152 carrés kilomètres.

À titre de vérification, les formules appliquées à toutes les cellules d'une grille de 720 par 360 couvrant la surface de la terre donnent une superficie totale de 4 * pi * (6371.0071809) ^ 2 kilomètres carrés, indiquant que le rayon authalique de la terre devrait être de 6371.0071809 kilomètres. Cela ne diffère de la valeur Wikipedia que dans le dernier chiffre significatif (environ un dixième de millimètre). (Je pense que les calculs de Wikipédia sont un peu décalés :-).

Comme vérifications supplémentaires, j'ai utilisé des versions de ces formules pour reproduire les annexes 4 et 5 dans Lev M. Bugayevskiy et John P. Snyder, Map Projections: A Reference Manual (Taylor & Francis, 1995). L'annexe 4 montre les longueurs d'arc de sections de 30 pieds de méridiens et de parallèles, exprimées au mètre près. Une vérification ponctuelle des résultats a montré un accord parfait. J'ai ensuite recréé la table avec des incréments de 0,0005 ', plutôt que des incréments de 0,5', et j'ai intégré numériquement les zones de quadrilatère estimées avec ces longueurs d'arc. La superficie totale de l'ellipsoïde a été reproduite avec précision à plus de huit chiffres significatifs. L'annexe 5 montre les valeurs de area(f)pour f = 0, 1/2, 1, ..., 90 degrés, multipliées par 1 / (2 * pi). Ces valeurs sont données au kilomètre carré le plus proche. Une vérification visuelle des valeurs proches de 0, 45 et 90 degrés a montré un accord parfait.

Cette formule exacte peut être appliquée en utilisant l'algèbre raster en commençant par une grille donnant les latitudes des limites supérieures de chaque cellule et une autre donnant les latitudes des limites inférieures. Chacun d'eux est, essentiellement, une grille de coordonnées y. (Dans chaque cas, vous souhaiterez peut-être créer sin(f)puis, puis zmet zpcomme résultats intermédiaires.) Soustraire les deux résultats, prendre la valeur absolue de cela et multiplier par la fraction q obtenue à la première étape (égale à 0,5 / 360 = 1/720 pour une largeur de cellule de 30 ', par exemple). Ce sera une grille dont les valeurs contiennent l' exactzones de chaque cellule (jusqu'à la précision numérique de la grille). Assurez-vous simplement d'exprimer les latitudes sous la forme attendue par la fonction sinus: de nombreuses calculatrices raster vous donneront des coordonnées en degrés mais attendez-vous aux radians pour leurs fonctions trigonométriques!

Pour mémoire, voici les zones exactes de 30 'par 30' cellules sur l'ellipsoïde WGS 84 de l'équateur jusqu'à un pôle, par intervalles de 30 ', à 11 chiffres (le même nombre utilisé pour le petit rayon 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

Les valeurs sont en kilomètres carrés.

Si vous souhaitez approximer ces zones ou simplement mieux comprendre leur comportement, la formule se réduit à une série de puissance suivant ce modèle:

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

où

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

(Une formule équivalente apparaît dans Bugayevskiy & Snyder, op. Cit. , Équation (2.1).)

Puisque e ^ 2 est si petit (environ 1/150 pour tous les modèles ellipsoïdaux de la terre) et z se situe entre 0 et 1, y est petit aussi. Ainsi, les termes y ^ 2, y ^ 3, ... deviennent rapidement plus petits, ajoutant plus de deux décimales de plus à chaque terme. Si nous devions ignorer tout à fait y , la formule serait celle de l'aire d'une sphère de rayon b . Les termes restants peuvent être compris comme corrigeant le renflement équatorial de la Terre.

Éditer

Certaines questions ont été soulevées concernant la façon dont un calcul de distance géodésique de surface se compare à ces formules exactes. La méthode de la distance géodésique rapproche chaque quadrilatère des géodésiques, plutôt que des parallèles, qui relient ses coins horizontalement, et applique la formule euclidienne pour un trapèze. Pour les petits quadrangles, tels que les quadrants de 30 ', ce paramètre est légèrement biaisé et a une précision relative comprise entre 6 et 10 parties par million. Voici un tracé de l'erreur pour WGS 84 (ou tout autre ellipsoïde terrestre raisonnable, d'ailleurs):

Ainsi, si (1) vous avez facilement accès aux calculs de distance géodésique et (2) peut tolérer une erreur de niveau ppm, vous pourriez envisager d'utiliser ces calculs géodésiques et de multiplier leurs résultats par 1,00000791 pour corriger le biais. Pour deux décimales de précision supplémentaires, soustrayez pi / 2 * cos (2f) / 10 ^ 6 du facteur de correction: le résultat sera précis à 0,04 ppm près.

La réponse à la question de radouxju dépend de la forme du pixel lorsqu'il est projeté sur l'ellipsoïde. Si le système de coordonnées du raster est la longitude et la latitude, alors le pixel est un rectangle de ligne de rhumb et la réponse de whuber peut être utilisée, ou, plus généralement, vous pouvez utiliser la formule pour un polygone dont les bords sont des lignes de rhumb. Si le système de coordonnées est une projection conforme à grande échelle (UTM, plan d'état, etc.), il serait plus précis d'approximer les arêtes par géodésiques et d'utiliser la formule pour un polygone géodésique. Les polygones géodésiques sont probablement les meilleurs pour une utilisation générale, car, contrairement aux polygones de ligne de rhumb, ils sont "bien comportés" près des pôles.

Les implémentations des formules pour les polygones géodésiques et rhumb sont fournies par ma bibliothèque GeographicLib . La zone géodésique est disponible en plusieurs langues; la zone de la ligne de rhumb est uniquement en C ++. Il existe une version en ligne (géodésique + rhumb line) disponible ici . La précision de ces calculs est généralement meilleure que 0,1 mètre carré.

Vous devrez juger sur crédible / officiel ... Les formules géodésiques sont dérivées de la zone sous la géodésique (Danielsen, 1989, abonnement requis), et des algorithmes pour la géodésique (Karney, 2013, libre accès). Les formules des lignes de rhumb sont données ici .

J'ai rencontré cette question en essayant de déterminer une formule pour la zone d'un pixel WGS84. Bien que la réponse de @ whuber contienne ces informations, il restait du travail à obtenir une formule pour la zone d'un pixel carré à une latitude donnée. J'ai inclus une fonction Python que j'ai écrite ci-dessous qui résume cela en un seul appel. Bien qu'il ne réponde pas directement à la question de l'affiche sur la zone d'un raster ENTIER (bien que l'on puisse additionner les zones de tous les pixels), je pense que ce sont des informations utiles pour quelqu'un qui recherche peut-être un calcul similaire.

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])