Cara yang lebih akurat untuk menghitung luas raster

Dalam pekerjaan sehari-hari saya, saya terus-menerus diminta untuk menghitung area dataset raster global dalam proyeksi geografis pada resolusi 30 arc detik. Kumpulan data ini biasanya merupakan hasil dari operasi Combine (contoh khasnya adalah kelas vegetasi yang dikombinasikan dengan layer country). Untuk melakukan ini, unit kami membuat dataset raster dengan luas setiap piksel dalam proyeksi geografis pada 30 arc detik. Dengan kisi area ini, zonalstat dilakukan untuk menjumlahkan area untuk setiap kelas. Karena saya tidak yakin bagaimana grid area ini dibuat, saya selalu bertanya-tanya apakah pendekatan ini lebih akurat hanya dengan memproyeksi ulang raster dalam proyeksi area yang sama (dari tes sederhana hasil dari kedua metode ini serupa). Apakah ada yang mengalami situasi yang sama?

Ada rumus tepat yang relatif sederhana untuk area segi empat bola yang dibatasi oleh paralel (garis lintang) dan garis meridian (garis bujur). Ia dapat diturunkan secara langsung menggunakan sifat-sifat dasar elips (sumbu utama a dan sumbu minor b ) yang diputar di sekitar sumbu minornya untuk menghasilkan ellipsoid. (Derivasi membuat latihan Kalkulus integral yang bagus tapi saya yakin akan sedikit menarik di situs ini.)

Rumus ini disederhanakan dengan memecah perhitungan menjadi langkah-langkah dasar.

Pertama, jarak antara batas timur dan barat - meridian l0 dan l1 - adalah sebagian kecil dari seluruh lingkaran yang sama dengan q = (l1 - l0) / 360 (ketika meridian diukur dalam derajat) atau 1 = ( l1 - l0) / (2 * pi) (ketika meridian diukur dalam radian). Temukan luas seluruh irisan yang terletak di antara paralel f0 dan f1 dan gandakan saja dengan q .

Kedua, kita akan menggunakan rumus untuk bidang irisan horizontal ellipsoid yang dibatasi oleh Khatulistiwa (pada f0 = 0) dan sejajar pada lintang f (= f1). Luas irisan antara dua garis lintang f0 dan f1 (terletak di belahan bumi yang sama) akan menjadi perbedaan antara area yang lebih besar dan lebih kecil.

Akhirnya, asalkan model tersebut benar-benar sebuah ellipsoid (dan bukan bola), area irisan antara Khatulistiwa dan paralel pada garis lintang f diberikan oleh

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

di mana adan bpanjang sumbu utama dan minor dari elips pembangkit, masing-masing,

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

adalah eksentrisitasnya, dan

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

(Ini jauh lebih sederhana daripada komputasi dengan geodesik, yang hanya merupakan perkiraan paralelnya saja. Harap perhatikan komentar oleh @cffk mengenai cara menghitung log(zp/zm)dengan cara yang menghindari hilangnya presisi pada lintang rendah.)

area(f) adalah area irisan buram dari garis khatulistiwa hingga lintang f (sekitar 30 derajat utara dalam ilustrasi. X dan Y adalah geosentris sumbu koordinat Kartesius yang ditunjukkan untuk referensi.

Untuk ellipsoid WGS 84, gunakan nilai-nilai konstan

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

menarik

e = 0.08181919084296

(Untuk model bola dengan a = b , rumus menjadi tidak pasti. Anda harus mengambil batas sebagai e -> 0 dari atas, yang kemudian direduksi menjadi rumus standar 2 * pi * a^2 * sin(f).)

Menurut rumus ini, segi empat 30 'oleh 30' berdasarkan Khatulistiwa memiliki luas 3077,2300079129 kilometer persegi, sedangkan segi empat 30 'oleh 30' menyentuh tiang (yang benar-benar hanya sebuah segitiga) memiliki luas hanya 13.6086152 persegi kilometer.

Sebagai tanda centang, rumus yang diterapkan pada semua sel kotak berukuran 720 x 360 yang menutupi permukaan bumi memberikan total luas permukaan 4 * pi * (6371.0071809) ^ 2 kilometer persegi, yang mengindikasikan bahwa jari - jari othalika bumi harus 6371.0071809 kilometer. Ini berbeda dari nilai Wikipedia hanya pada angka signifikan terakhir (sekitar sepersepuluh milimeter). (Saya pikir perhitungan Wikipedia sedikit tidak aktif :-).

Sebagai pemeriksaan tambahan, saya menggunakan versi formula ini untuk mereproduksi Lampiran 4 dan 5 di Lev M. Bugayevskiy & John P. Snyder, Proyeksi Peta: Manual Referensi (Taylor & Francis, 1995). Lampiran 4 menunjukkan panjang busur bagian meridian dan paralel sepanjang 30'-panjang, diberikan ke meter terdekat. Hasil pemeriksaan langsung menunjukkan persetujuan yang sempurna. Saya kemudian menciptakan kembali tabel dengan kenaikan 0,0005 ', bukan kenaikan 0,5', dan secara numerik mengintegrasikan bidang segi empat seperti yang diperkirakan dengan panjang busur ini. Total area ellipsoid direproduksi secara akurat hingga lebih baik dari delapan angka penting. Lampiran 5 menunjukkan nilai area(f)untuk f = 0, 1/2, 1, ..., 90 derajat, dikalikan dengan 1 / (2 * pi). Nilai-nilai ini diberikan kepada kilometer persegi terdekat. Pemeriksaan visual dari nilai di dekat 0, 45 dan 90 derajat menunjukkan persetujuan sempurna.

Rumus yang tepat ini dapat diterapkan dengan menggunakan aljabar raster yang dimulai dengan kisi yang memberi garis lintang batas atas setiap sel dan sel lainnya memberi garis lintang batas bawah. Masing-masing, pada dasarnya, adalah kisi-kisi y. (Dalam setiap kasus Anda mungkin ingin membuat sin(f)dan kemudian zmdan zpsebagai hasil antara.) Kurangi dua hasil, ambil nilai absolut dari itu, dan kalikan dengan fraksi q yang diperoleh pada langkah pertama (sama dengan 0,5 / 360 = 1/720 untuk lebar sel 30 ', misalnya). Ini akan menjadi kotak yang nilainya berisi angka yang tepatarea masing-masing sel (hingga presisi numerik grid sendiri). Pastikan untuk mengekspresikan garis lintang dalam bentuk yang diharapkan oleh fungsi sinus: banyak kalkulator raster akan memberi Anda koordinat dalam derajat tetapi mengharapkan radian untuk fungsi trigonometri mereka!

Sebagai catatan, berikut adalah area tepat sel 30 'oleh 30' pada ellipsoid WGS 84 dari Ekuator hingga sebuah kutub, dalam interval 30 ', hingga 11 angka (angka yang sama digunakan untuk jari-jari kecil 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

Nilainya dalam kilometer persegi.

Jika Anda ingin memperkirakan area-area ini atau hanya memahami perilakunya dengan lebih baik, rumusnya dikurangi menjadi serangkaian daya yang mengikuti pola ini:

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

dimana

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

(Rumus yang setara muncul di Bugayevskiy & Snyder, op. Cit. , Persamaan (2.1).)

Karena e ^ 2 sangat kecil (sekitar 1/150 untuk semua model ellipsoidal di bumi) dan z terletak di antara 0 dan 1, y juga kecil. Jadi istilah y ^ 2, y ^ 3, ... dengan cepat menjadi lebih kecil, menambahkan lebih dari dua tempat desimal lebih presisi dengan setiap istilah. Jika kita mengabaikan y sama sekali, rumusnya adalah luas bidang radius b . Istilah yang tersisa dapat dipahami sebagai mengoreksi tonjolan khatulistiwa bumi.

Edit

Beberapa pertanyaan telah diajukan mengenai bagaimana perhitungan jarak geodesi pada area dibandingkan dengan formula yang tepat ini. Metode jarak geodesik mendekati setiap segi empat dengan geodesik, bukan paralel, yang menghubungkan sudut-sudutnya secara horizontal, dan menerapkan rumus Euclidean untuk trapesium. Untuk quadrangles kecil, seperti quads 30 ', ini bias sedikit rendah dan memiliki akurasi relatif antara 6 dan 10 bagian per juta. Berikut adalah plot kesalahan untuk WGS 84 (atau ellipsoid bumi yang masuk akal, dalam hal ini):

Jadi, jika (1) Anda memiliki akses mudah ke perhitungan jarak geodesik dan (2) dapat mentolerir kesalahan level ppm, Anda dapat mempertimbangkan menggunakan perhitungan geodesik tersebut dan mengalikan hasilnya dengan 1,00000791 untuk mengoreksi bias. Untuk dua tempat desimal yang lebih presisi, kurangi pi / 2 * cos (2f) / 10 ^ 6 dari faktor koreksi: hasilnya akan akurat hingga 0,04 ppm.

Jawaban untuk pertanyaan radouxju tergantung pada bentuk piksel ketika diproyeksikan ke ellipsoid. Jika sistem koordinat raster adalah bujur dan lintang, maka pikselnya adalah persegi panjang garis rhumb dan jawaban whuber dapat digunakan, atau, secara lebih umum, Anda dapat menggunakan rumus untuk poligon yang tepinya adalah garis rhumb. Jika sistem koordinat adalah proyeksi konformal skala besar (UTM, state plane, dll.), Akan lebih akurat untuk memperkirakan tepi dengan geodesik dan menggunakan rumus untuk poligon geodesik. Poligon geodesik mungkin paling baik untuk penggunaan umum, karena, tidak seperti poligon rhumb, poligon "berperilaku baik" dekat dengan kutub.

Implementasi formula untuk poligon garis geodesik dan rhumb disediakan oleh perpustakaan saya GeographicLib . Area geodesik tersedia dalam beberapa bahasa; area baris rhumb adalah C ++ saja. Ada versi online (geodesik + baris rhumb) tersedia di sini . Keakuratan perhitungan ini biasanya lebih baik dari 0,1 meter persegi.

Anda harus menilai berdasarkan kredibilitas / resmi ... Rumus geodesik diturunkan di Area di bawah geodesik (Danielsen, 1989, diperlukan berlangganan), dan Algoritma untuk geodesik (Karney, 2013, akses terbuka). Rumus baris rhumb diberikan di sini .

Saya menemukan pertanyaan ini ketika mencoba menentukan formula untuk area piksel WGS84. Sementara jawaban @ whuber memang mengandung informasi ini, masih ada beberapa pekerjaan untuk mendapatkan formula untuk area piksel derajat persegi pada garis lintang tertentu. Saya telah menyertakan fungsi Python yang saya tulis di bawah ini yang mengabstraksi ini menjadi satu panggilan. Meskipun tidak secara langsung menjawab pertanyaan poster tentang area seluruh raster (walaupun orang dapat menjumlahkan area semua piksel), saya pikir itu masih informasi yang berguna bagi seseorang yang mungkin mencari perhitungan yang sama.

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