Membuat DEM skala nano dengan GDAL

14

Mungkin sedikit pertanyaan aneh, tapi izinkan saya memberi Anda penjelasan singkat tentang latar belakang sebelum pertanyaan aktual saya:

Atomic force microscopy (AFM) adalah metode yang, singkatnya (dan setahu saya) memungkinkan peneliti memindai area di skala mikro dan nano. Ini bekerja dengan "memindai" suatu area menggunakan semacam probe. Lebih banyak yang sulit untuk saya jelaskan, karena saya tidak memiliki pemahaman yang nyata tentang hal itu. Apa yang saya tahu, dan apa yang memicu keingintahuan saya adalah bahwa hasilnya sebenarnya adalah sebuah "kotak" nilai "tinggi" (sebuah matriks katakanlah nilai 512x512 yang menggambarkan ketinggian penyelidikan pada saat itu).

Saya kemudian berpikir: Yah, terlepas dari skala, ini sebenarnya adalah model elevasi digital! Dan, ini berarti bahwa jika saya dapat mengatur untuk membuat file DEM sebagaimana dipahami oleh alat GIS saya dapat menerapkan analisis GIS untuk itu!

Karena itu, pekerjaan penting saya yang lain di lab yang memiliki mesin AFM, dan menggunakannya di salah satu proyeknya. Saya mendapatkan beberapa file pindaian darinya, dan telah berhasil, menggunakan Python (struct dan numpy), untuk mengurai file-file biner ini dan apa yang saya miliki sekarang adalah array numpy dengan ukuran 512x512 yang diisi dengan nilai-nilai int16.

Apa yang saya rencanakan selanjutnya, dan apa yang saya perlu bantuan adalah bagian "pemetaan ke DEM yang tepat". Saya memiliki pengetahuan tentang DEMS, tetapi ketika sampai pada generasi mereka yang sebenarnya, saya cukup baru.

Yang saya pikirkan adalah saya harus melakukan georeferensi data saya dengan cara tertentu, dan untuk ini saya memerlukan sistem koordinat custom (planar). Saya membayangkan sistem koordinat saya akan menggunakan satuan mikro atau nano meter. Maka itu hanya masalah menemukan ukuran area yang dipindai dengan AFM (ini saya percaya ada di suatu tempat dalam file biner, anggap ini diketahui).

pembaruan : Saya juga memiliki beberapa pemindaian pada resolusi yang berbeda, tetapi pada area yang sama. Sebagai contoh, saya memiliki informasi ini tentang dua pemindaian:

gambar yang lebih besar:

Scan Size: 51443.5 nm
X Offset: 0 nm
Y Offset: 0 nm

gambar lebih kecil (detail):

Scan Size: 5907.44 nm
X Offset: 8776.47 nm
Y Offset: 1486.78 nm

Apa yang saya pikirkan adalah bahwa sistem koordinat kustom saya harus berasal dari 0,0 dan untuk gambar yang lebih besar saya dapat menetapkan piksel 0,0 nilai koordinat (0,0) dan piksel 512,512 nilai koordinat (51443,5, 51443,5 ) (Kira Anda mendapatkan gambar untuk poin lain yang dibutuhkan).

Kemudian, gambar yang lebih besar akan memetakan piksel (0,0) hingga (8776,47, 1486,78) dan (512,512) hingga (8776,47 + 5907,44, 1486,78 + 5907,44)

Pertanyaan 1 adalah : Bagaimana cara membuat proj4 def untuk sistem koordinat seperti itu? Yaitu: bagaimana cara menetapkan "koordinat dunia nyata" ini ke sistem koordinat kustom saya (atau, jika saya mengikuti saran pengacak dan menggunakan sistem koordinat lokal dan berbohong tentang unit (yaitu memperlakukan nanometer saya sebagai kilometer)

Kemudian saya harus mentransfer Array 2 dimensi saya yang numpy ke format file DEM Georeferensi. Saya berpikir untuk menggunakan GDAL (atau, lebih tepatnya pengikatan Python).

Pertanyaan ke-2 adalah : Bagaimana cara membuat DEM georeferensi dari data "sewenang-wenang" seperti milik saya? Lebih disukai dalam Python dan menggunakan perpustakaan open source.

Sisanya kemudian harus cukup mudah, hanya masalah menggunakan alat analisis yang tepat. Masalahnya adalah, tugas ini didorong oleh keingintahuan saya sendiri, jadi saya tidak yakin apa yang harus saya lakukan dengan DEM skala nano. Ini memohon

Pertanyaan 3 : Apa yang harus dilakukan dengan DEM skala nano? Analisis macam apa yang dapat dilakukan, alat apa yang tepat untuk analisis DEM dan akhirnya: apakah layak membuat peta dengan garis bukit dan garis kontur dari data ini? :)

Saya menyambut semua saran dan petunjuk, tetapi perlu diingat bahwa saya sedang mencari alternatif gratis, karena ini adalah proyek yang benar-benar berbasis hobi tanpa anggaran atau dana (dan saya tidak memiliki akses ke aplikasi GIS yang dilisensikan). Selain itu, saya tahu bahwa Bruker, perusahaan yang menjual mesin AFM ini, mengirimkan beberapa perangkat lunak, tetapi menggunakan itu tidak akan menyenangkan.

atlefren
sumber
Menyenangkan dan menarik! Bisakah Anda memposting beberapa data sampel? Apakah itu persyaratan untuk memiliki skala nanometer pada proyeksi? Berpikir mungkin lebih mudah untuk ditingkatkan, meskipun itu "curang" sedikit. Btw saya kira Anda bisa datang jauh dengan GDAL / ogr, meskipun masalah proyeksi masih perlu ditangani. gdal.org/gdal_grid.html
alexanno
Terima kasih! Kira ini lebih merupakan komentar daripada jawaban. Sejauh skala nanometer saya akan mengatakan bahwa apa pun yang bekerja sangat bagus, tetapi peoection skala nano yang benar akan menjadi yang paling keren. Ketika datang ke sampel data saya harus memeriksa apakah ada beberapa batasan, tetapi pada dasarnya ini adalah 2 dim matrix nilai int16.
atlefren
3
Mengapa Anda membutuhkan parameter proj4? Anda tidak akan mentransfer data ini ke sistem koordinat lain (khususnya geografis). Imho Anda harus dapat melakukan semua analisis Anda tanpa sistem koordinat. Apa yang kamu mengerti dengan DEM? Ada beberapa jenis seperti permukaan triangulasi (mis. Lakukan delauny triangulasi) atau peta raster (Anda sudah memilikinya). Ini tentu saja sangat tergantung pada perangkat lunak analisis. Tentu saja Anda dapat membuat peta lain jika diperlukan sebagai output untuk memahami probe. Anda mungkin melihat code.google.com/p/mtex untuk analisis gandum.
mistapink
3
Alasan saya (berpikir) saya membutuhkan CRS adalah ini: jika saya hanya membuat, katakanlah GeoTIFF tanpa menetapkan CRS, satuan ukurannya adalah piksel. Bagaimana jika saya ingin mengukur jarak? Dan, bagaimana jika saya memiliki dua pemindaian AFM dan tahu bagaimana mereka saling berhubungan (dalam hal skala dan mengimbangi dari beberapa titik).
Menugaskan
1
Saya biasanya mengatasi data tersebut dengan mengatur sistem koordinat lokal (dengan asal yang bertepatan dengan asal gambar) dan "berbohong" tentang unit. Misalnya, Anda dapat menetapkan bahwa unit adalah kilometer ketika mereka benar-benar nanometer, yang membuatnya mudah secara mental untuk mengkonversi bolak-balik. Tentu saja Anda tidak akan melakukan proyeksi ulang, jadi itu bukan masalah. Menyiapkan sistem koordinat ini sama dengan melakukan georeferensi setiap DEM; itu bisa sesederhana membuat file dunia yang tidak diputar.
whuber

Jawaban:

4

Yah, sepertinya saya menyelesaikan masalah 1 & 2 setidaknya. Kode sumber lengkap di github , tetapi beberapa penjelasan di sini:

Untuk CRS khusus, saya memutuskan (sesuai saran Whubers) untuk "menipu" dan menggunakan meter sebagai unit. Saya menemukan "cr lokal" di apatialreference.org ( SR-ORG: 6707 ):

LOCAL_CS["Non-Earth (Meter)",
    LOCAL_DATUM["Local Datum",0],
    UNIT["Meter",1.0],
    AXIS["X",NORTH],
    AXIS["Y",EAST]]

Menggunakan Python dan GDAL ini cukup mudah dibaca:

def get_coordsys():
    #load our custom crs
    prj_text = open("coordsys.wkt", 'r').read()
    srs = osr.SpatialReference()
    if srs.ImportFromWkt(prj_text):
        raise ValueError("Error importing PRJ information" )

    return srs

Juga, memanaskan DEM dengan GDAL sebenarnya cukup mudah (saya berakhir dengan geo-tiff single-band). Baris parser.read_layer (0) mengembalikan matriks 512x512 yang saya jelaskan sebelumnya.

def create_dem(afmfile, outfile):

    #parse the afm data
    parser = AFMParser(afmfile)

    #flip to account for the fact that the matrix is top-left, but gdal is bottom-left
    data = flipud(parser.read_layer(0))

    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(
        outfile,
        data.shape[1],
        data.shape[0],
        1 ,
        gdal.GDT_Int16 ,
    )

    dst_ds.SetGeoTransform(get_transform(parser, data))
    dst_ds.SetProjection(get_coordsys().ExportToWkt())
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds = None

Bagian triockiest mencari tahu cara "georeferensi" file saya dengan benar, saya akhirnya menggunakan SetGeoTransform , mendapatkan parameter sebagai berikut:

def get_transform(parser, data):
    #calculate dims
    scan_size, x_offset, y_offset = parser.get_size()
    x_size = data.shape[0]
    y_size = data.shape[1]
    x_res = scan_size / x_size
    y_res = scan_size / y_size

    #set the transform (offsets doesn't seem to work the way I think)
    #top left x, w-e pixel resolution, rotation, 0 if image is "north up", top left y, rotation, 0 if image is "north up", n-s pixel resolution
    return [0, x_res, 0, 0, 0, y_res]

Bagian terakhir ini mungkin adalah bagian yang paling tidak saya yakini, yang sebenarnya saya cari adalah baris * gdal_transform -ullr *, tetapi saya tidak dapat menemukan cara untuk melakukan ini secara terprogram.

Saya dapat membuka GeoTIFF saya di Qgis dan melihatnya (dan membandingkannya secara visual dengan hasil dari program Bruker yang terlihat benar), tetapi saya belum benar-benar menjawab pertanyaan saya 3; apa yang harus dilakukan dengan data ini. Jadi, di sini saya terbuka untuk saran!

atlefren
sumber
Salah satu tantangan yang menarik adalah membandingkan jarak pada DEM dengan jarak antara tempat-tempat di dunia untuk memberi pemirsa gambaran betapa kecilnya skala nano. Misalnya htwins.net/scale2
blah238