Przycinanie rastra warstwą wektorową za pomocą GDAL

Zainstalowałem GDAL za pomocą instalatora Osgeo. Jak mogę programowo przyciąć warstwę rastrową warstwą wektorową? Czy jest dostępny interfejs API GDAL, który może mi w tym pomóc? Używam Pythona.

Nie jestem pewien co do interfejsu gdal, jest void* GDALWarpOptions::hCutlinew Warp Options, do którego odwołuje się samouczek API Warp API , ale nie ma wyraźnych przykładów. Czy na pewno potrzebujesz programowej odpowiedzi? Narzędzia wiersza poleceń mogą to zrobić natychmiast po wyjęciu z pudełka:

1. utwórz plik kształtu zawierający tylko obszar zainteresowania przycinający wielokąt
2. służy ogrinfodo określania zasięgu przycinającego pliku kształtu
3. służy gdal_translatedo przycinania do zakresu kształtu
4. użyj gdalwarpz -cutlineparametrem

Kroki 2 i 3 są do optymalizacji, z którą można sobie poradzić gdalwarp -cutline ....

Zobacz Przycinanie rastrów za pomocą GDAL przy użyciu wielokątów z Linfinity dla rozwiązań opartych na systemie Linux, wszystkie opakowane w jeden skrypt. Inny przykładowy przykład można znaleźć w samouczku Michaela Coreya dotyczącym tworzenia cieni dla Mapnika .

Joel Lawhead z GeospatialPython ma kompletny przykład python w Clip rastrze, używając shapefile , dobrze napisanego samouczka. Konieczne będzie zainstalowanie biblioteki obrazów Python (PIL), która nie jest zawarta w Osgeo4W (dla której może być konieczne dodanie Pythona o4w do rejestru systemu Windows, aby program instalacyjny działał).

Wydaje się, że temat ten zawsze powraca. Sam nie wiedziałem, że GDAL> 1.8 jest tak zaawansowany, że już zapewnia ci uczciwą obsługę wiersza poleceń, aby wykonać to zadanie.

Komentarz Mike'a Toews jest całkiem przydatny, ale możesz po prostu zrobić na przykład:

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

Możesz zawinąć to polecenie w skrypt Pythona za pomocą doskonałego modułu podprocesu .

Jedną z rzeczy, która była dla mnie naprawdę problematyczna, było to, że musiałem dostarczyć minimalne rozwiązanie tego problemu, co oznacza, że ​​jest tak proste, jak to możliwe i nie wymaga wielu zewnętrznych zależności. Korzystanie z Python Imaging Library jak w samouczku Joela Lawhead'a jest fajne, ale wpadłem na następujące rozwiązanie: używając tablic zamaskowanych Numpy.
Nie wiem, czy to jest lepsze, ale to wiedziałem wtedy (3 lata temu ...).
Pierwotnie utworzyłem prawidłowy obszar danych wewnątrz oryginalnego rastra (np. Zakres rastra wyjściowego, gdzie to samo), ale podobał mi się pomysł zmniejszenia rastra również (np. -Crop_to_cutline), więc zacząłem world2Pixelod Joela Lawhead. Oto moje własne rozwiązanie:

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)

pełny opis class MaskRasteri jego metod można znaleźć na githubie mojego projektu .

Używając tego kodu nadal będziesz musiał używać GDAL. Jednak planuję używać w przyszłości czystego Pythona tam, gdzie mogę, ponieważ docelowi odbiorcy mojego oprogramowania mają trudności ze zbyt wieloma zależnościami (używam Debiana do tworzenia oprogramowania, a klienci używają Windows 7 ...).