Przycinanie rastra warstwą wektorową za pomocą GDAL

26

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.

Vicky
źródło

Odpowiedzi:

13

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 .

matowe wilkie
źródło
Matt, możesz sobie przypomnieć trac.osgeo.org/gdal/ticket/1599 wygląda na to, że cutline to spełnia
Mike T
10

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

Oz123
źródło
Podoba mi się podany przykład wiersza polecenia, ale czy możesz wyjaśnić, co robi argument -crop_to_cutline? Nie jestem pewien, jaki jest cel przycinania pliku kształtów określany przez -cutline.
hendra
1
opcja -cutline przycina raster do wewnętrznej ramki granicznej warstwy wielokąta. Np. Jeśli jest mniejszy w zakresie, raster wyjściowy również będzie mniejszy. Bez tego raster wyjściowy będzie miał taki sam rozmiar jak oryginał, ale z NULL we wszystkich punktach poza obszarem zainteresowania.
Oz123