Kliping raster dengan layer vektor menggunakan GDAL

Saya telah menginstal GDAL menggunakan installer Osgeo. Bagaimana saya bisa memotong layer raster dengan layer vektor secara terprogram? Apakah ada API GDAL yang tersedia yang dapat membantu saya dalam hal ini? Saya menggunakan Python.

Saya tidak yakin tentang api gdal, ada void* GDALWarpOptions::hCutlinedi Opsi Warp direferensikan dari tutorial Warp API , tetapi tidak ada contoh eksplisit. Apakah Anda yakin Anda membutuhkan jawaban terprogram? Utilitas baris perintah dapat melakukannya di luar kotak:

1. buat shapefile yang hanya berisi bidang kliping poligon yang menarik
2. gunakan ogrinfountuk menentukan sejauh mana shapefile kliping
3. gunakan gdal_translateuntuk klip ke bentuk luasan
4. gunakan gdalwarpdengan -cutlineparameter

Langkah 2 & 3 adalah untuk optimasi, Anda bisa bertahan dengan adil gdalwarp -cutline ....

Lihat Memotong raster dengan GDAL menggunakan poligon dari Linfinity untuk solusi berbasis linux semua dibungkus dalam satu skrip. Contoh cutline lain dapat dilihat pada tutorial Michael Corey yang menciptakan hillshades untuk Mapnik .

Joel Lawhead dari GeospatialPython memiliki contoh python lengkap di Clip raster menggunakan shapefile , tutorial yang ditulis dengan baik. Anda harus menginstal Python Image Library (PIL) yang tidak termasuk dalam Osgeo4W (untuk itu Anda mungkin perlu menambahkan o4w python ke registri windows untuk mendapatkan program instal agar berfungsi).

Tampaknya subjek ini selalu kembali. Saya sendiri tidak tahu bahwa GDAL> 1.8 sangat canggih sehingga sudah memberi Anda penanganan baris perintah yang adil untuk melakukan tugas itu.

Komentar dari Mike Toews cukup berguna tetapi Anda bisa melakukannya misalnya:

gdalwarp -of GTiff -cutline DATA/area_of_interest.shp -cl area_of_interest  -crop_to_cutline DATA/PCE_in_gw.asc  data_masked7.tiff 

Anda bisa membungkus perintah ini di dalam skrip python dengan modul subproses yang sangat baik .

Satu hal yang benar-benar bermasalah bagi saya adalah bahwa saya perlu menyediakan solusi minimal untuk masalah itu, yang berarti sesederhana mungkin dan tidak memerlukan banyak ketergantungan eksternal. Penggunaan Python Imaging Library seperti dalam tutorial oleh Joel Lawhead rapi, tapi saya datang dengan solusi berikut: menggunakan array bertopeng Numpy.
Saya tidak tahu apakah itu lebih baik, tapi itu yang saya tahu (3 tahun lalu ...).
Awalnya saya membuat area data yang valid di dalam raster asli (misalnya sejauh mana output raster sama), tetapi saya menyukai gagasan membuat raster juga lebih kecil (misalnya -crop_to_cutline), jadi saya mengadopsi world2Pixeldari Joel Lawhead. Ini solusi saya sendiri:

def RasterClipper():
    craster = MaskRaster()
    contraster2 = 'PCE_in_gw.aux'
    craster.reader("DATA/"+contraster2.replace('aux','asc'))
    xres, yres = craster.extent[1], craster.extent[1]
    craster.fillrasterpoints(xres, yres)
    craster.getareaofinterest("DATA/area_of_interest.shp")
    minX, maxX=craster.new_extent [0]-5,craster.new_extent[1]+5
    minY, maxY= craster.new_extent [2]-5,craster.new_extent[3]+5
    ulX, ulY=world2Pixel(craster.extent, minX, maxY)
    lrX, lrY=world2Pixel(craster.extent, maxX, minY)
    craster.getmask(craster.corners)
    craster.mask=np.logical_not(craster.mask)
    craster.mask.resize(craster.Yrange.size,craster.Xrange.size)
    # choose all data points inside the square boundaries of the AOI,
    # replace all other points with NULL
    craster.cdata= np.choose(np.flipud(craster.mask), (craster.data, -9999))
    # resise the data set to be the size of the squared polygon
    craster.ccdata=craster.cdata[ulY:lrY, ulX:lrX]
    craster.writer("ccdata2m.asc",craster.ccdata, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)
    # in second step we rechoose all the data points which are inside the
    # bounding vertices of AOI
    # need to re-define our raster points
    craster.xllcorner, craster.yllcorner = minX, minY
    craster.xurcorner, craster.yurcorner = maxX, maxY
    craster.fillrasterpoints(10,10)
    craster.getmask(craster.boundingvertices) # just a wrapper around matplotlib.nxutils.points_in_poly
    craster.data=craster.ccdata
    craster.clip2(new_extent_polygon=craster.boundingvertices)
    craster.data = np.ma.MaskedArray(craster.data, mask=craster.mask)
    craster.data = np.ma.filled(craster.data, fill_value=-9999)
    # write the raster to disk
    craster.writer("ccdata2m_clipped.asc",craster.data, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)

untuk penjelasan lengkap tentang class MaskRasterdan metodenya, lihat github proyek saya .

Menggunakan kode ini Anda masih harus menggunakan GDAL. Namun, rencananya adalah untuk menggunakan Python murni di masa depan di mana saya bisa, karena audiens yang dituju perangkat lunak saya mengalami kesulitan dengan terlalu banyak dependensi (saya menggunakan Debian untuk mengembangkan perangkat lunak, dan klien menggunakan Windows 7 ...).