Menemukan lokasi dengan nilai tertinggi dalam raster menggunakan ArcGIS Desktop?

12

Menggunakan ArcGIS 10, saya memiliki raster di mana saya ingin mencari pixel dengan nilai maksimum di raster dan mengembalikan lokasinya (pusat pixel) dalam derajat desimal. Saya ingin mengulangi melalui proses ini mengembalikan lokasi nilai tertinggi kedua dari raster, kemudian ketiga, dan seterusnya sehingga pada akhirnya saya memiliki daftar lokasi N yang memiliki nilai tertinggi dalam raster secara berurutan.

Saya membayangkan bahwa ini mungkin paling mudah dilakukan dengan menggunakan skrip Python tapi saya terbuka untuk ide lain jika ada cara yang lebih baik.

mga
sumber
sudahkah Anda mencoba mengonversi kisi ke poin kemudian menambahkan bidang X, Y dan menyortir?
Jakub Sisak GeoGraphics
Apakah nilai raster mengapung atau bilangan bulat?
whuber
@ Jakub - Tidak, belum. Saya mungkin hanya akan tertarik pada 1% atau lebih dari poin jadi saya tidak tahu apakah perlu menambahkan x, y bidang untuk semua poin kemudian mengurutkan. Mungkin jika tidak ada pilihan yang lebih efisien?
mga
@whuber - Nilai raster adalah float.
mga
@mga, patut dicoba. Konversi cukup cepat dan menambahkan XY juga merupakan alat default. Menghapus catatan yang tidak diinginkan adalah operasi lurus ke depan dan semua bisa digabungkan menjadi satu model. Hanya sebuah ide.
Jakub Sisak GeoGraphics

Jawaban:

5

Jika Anda senang menggunakan R , ada paket yang disebut raster . Anda dapat membaca dalam raster menggunakan perintah berikut:

install.packages('raster')
library(raster)
test <- raster('F:/myraster')

Kemudian, saat Anda melihatnya (dengan mengetik test), Anda dapat melihat info berikut:

class       : RasterLayer 
dimensions  : 494, 427, 210938  (nrow, ncol, ncell)
resolution  : 200, 200  (x, y)
extent      : 1022155, 1107555, 1220237, 1319037  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
values      : F:/myraster 
min value   : 0 
max value   : 1 

Mungkin ada cara yang lebih baik untuk memanipulasi raster, tetapi salah satu cara menemukan info yang Anda inginkan adalah menemukan nilai tertinggi, dan mendapatkan lokasi matriksnya, lalu menambahkannya ke tingkat yang lebih rendah.

djq
sumber
1
+1 Untuk referensi ini. Setelah membaca raster R, Anda dapat menggunakan Rfungsi standar atau getValuesmetode untuk mengakses nilai sel. Dari sana mudah untuk mengidentifikasi nilai tertinggi dan lokasi mereka.
whuber
1
Berkat rekomendasi Anda, inilah yang akhirnya saya lakukan. Menggunakan paket raster di R sangat mudah dibandingkan dengan mencoba ini di ArcGIS. Saya juga terus menggunakan analisis spasial lainnya di R dan sangat senang dengan hasilnya. Saran bagus!
mga
8

Jawabannya dapat diperoleh dengan menggabungkan kisi indikator dari 1% nilai teratas dengan kisi lintang dan bujur. Kuncinya terletak pada pembuatan kisi indikator ini, karena ArcGIS (masih! Setelah 40 tahun!) Tidak memiliki prosedur untuk memeringkat data raster.

Salah satu solusi untuk raster floating-point adalah iteratif, tetapi untungnya cepat . Biarkan n menjadi jumlah sel data. The distribusi kumulatif empiris dari nilai-nilai terdiri dari semua pasangan (z, n (z)) di mana z adalah nilai dalam grid dan n (z) adalah jumlah sel di grid dengan nilai kurang dari atau sama dengan z . Kita mendapatkan kurva yang menghubungkan (-infinity, 0) ke (+ infinity, n) dari urutan simpul-simpul ini yang dipesan oleh z . Dengan demikian mendefinisikan fungsi f , di mana (z, f (z)) selalu terletak pada kurva. Anda ingin menemukan titik (z0, 0,99 * n) pada kurva ini.

Dengan kata lain, tugasnya adalah menemukan nol dari f (z) - (1-0.01) * n . Lakukan ini dengan rutinitas nol-temuan (yang dapat menangani fungsi sewenang-wenang: yang ini tidak dapat dibedakan). Yang paling sederhana, yang seringkali efisien, adalah menebak-dan-periksa: awalnya Anda tahu z0 terletak antara nilai minimum zMin dan zMax maksimum. Tebak nilai wajar apa pun di antara keduanya. Jika tebakannya terlalu rendah, atur zMin = z0; jika tidak, atur zMax = z0. Sekarang ulangi. Anda akan segera menyatu dengan solusi; Anda cukup dekat ketika zMax dan zMin cukup dekat. Agar konservatif, pilih nilai akhir zMin sebagai solusinya: mungkin mengambil beberapa poin tambahan yang dapat Anda buang nanti. Untuk pendekatan yang lebih canggih, lihat Bab 9 dari Resep Numerik (tautan menuju ke versi gratis yang lebih lama).

Melihat kembali algoritma ini menunjukkan Anda hanya perlu melakukan dua jenis operasi raster : (1) pilih semua sel kurang dari atau sama dengan beberapa nilai target dan (2) hitung sel yang dipilih. Itu adalah salah satu operasi paling sederhana dan tercepat di sekitar. (2) dapat diperoleh sebagai hitungan zonal atau dengan membaca satu catatan dari tabel atribut dari kisi pilihan.

whuber
sumber
7

Saya melakukan ini beberapa waktu yang lalu, meskipun solusi saya menggunakan GDAL (jadi, ini bukan hanya untuk ArcGIS). Saya pikir Anda bisa mendapatkan array NumPy dari raster di ArcGIS 10, tapi saya tidak tahu pasti. NumPy memasok pengindeksan array yang sederhana dan kuat, seperti argsortdan lainnya. Contoh ini tidak menangani NODATA atau mengubah koordinat dari yang diproyeksikan ke lat / long (tapi ini tidak sulit dilakukan dengan osgeo.osr, disediakan dengan GDAL)

import numpy as np
from osgeo import gdal

# Open raster file, and get GeoTransform
rast_src = gdal.Open(rast_fname)
rast_gt = rast_src.GetGeoTransform()

def get_xy(r, c):
    '''Get (x, y) raster centre coordinate at row, column'''
    x0, dx, rx, y0, ry, dy = rast_gt
    return(x0 + r*dx + dx/2.0, y0 + c*dy + dy/2.0)

# Get first raster band
rast_band = rast_src.GetRasterBand(1)

# Retrieve as NumPy array to do the serious work
rast = rast_band.ReadAsArray()

# Sort raster pixels from highest to lowest
sorted_ind = rast.argsort(axis=None)[::-1]

# Show highest top 10 values
for ind in sorted_ind[:10]:
    # Get row, column for index
    r, c = np.unravel_index(ind, rast.shape)
    # Get [projected] X and Y coordinates
    x, y = get_xy(r, c)
    print('[%3i, %3i] (%.3f, %.3f) = %.3f'%
          (r, c, x, y, rast[r, c]))

Menampilkan berikut ini untuk file raster pengujian saya:

[467, 169] (2813700.000, 6353100.000) = 844.538
[467, 168] (2813700.000, 6353200.000) = 841.067
[469, 168] (2813900.000, 6353200.000) = 840.705
[468, 168] (2813800.000, 6353200.000) = 840.192
[470, 167] (2814000.000, 6353300.000) = 837.063
[468, 169] (2813800.000, 6353100.000) = 837.063
[482, 166] (2815200.000, 6353400.000) = 833.038
[469, 167] (2813900.000, 6353300.000) = 832.825
[451, 181] (2812100.000, 6351900.000) = 828.064
[469, 169] (2813900.000, 6353100.000) = 827.514
Mike T
sumber
+1 Terima kasih telah berbagi ini. Saya tidak melihat ketidakmampuan untuk menangani NoData sebagai batasan: cukup konversi semua NoData ke nilai yang sangat negatif sebelum melanjutkan. Perhatikan juga, bahwa jika ada proyeksi ulang kisi - kisi terjadi, jawabannya kemungkinan akan berubah karena pemasangan ulang kisi-kisi, jadi biasanya orang tidak ingin proyeksi ulang terjadi secara otomatis selama perhitungan seperti ini. Sebagai gantinya, koordinat yang dilaporkan dapat diproyeksikan ulang sesudahnya. Dengan demikian solusi Anda sangat umum.
whuber
Penanganan NODATA dapat diimplementasikan dengan terlebih dahulu mendapatkan nilai dari raster NODATA = rast_band.GetNoDataValue(), kemudian dengan menggunakan nilai NaN ( rast[rast == NODATA] = np.nan) atau dengan menggunakan array masked ( rast = np.ma.array(rast, mask=(rast == NODATA))). Trik yang lebih rumit adalah argsortentah bagaimana menghapus nilai-nilai NODATA dari analisis, atau cukup lewati saja di for-loop jika mereka NaN / masked.
Mike T