Cortando raster com camada vetorial usando GDAL

26

Instalei o GDAL usando o instalador do Osgeo. Como cortar uma camada raster com uma camada vetorial programaticamente? Existe alguma API GDAL disponível que possa me ajudar com isso? Eu estou usando Python.

Vicky
fonte

Respostas:

13

Não tenho certeza sobre o gdal api, há void* GDALWarpOptions::hCutlinenas Opções de distorção referenciadas no tutorial da API de distorção , mas não há exemplos explícitos. Tem certeza de que precisa de uma resposta programática? Os utilitários de linha de comando podem fazer isso imediatamente:

  1. crie um shapefile contendo apenas a área de interesse do polígono de recorte
  2. use ogrinfopara determinar a extensão do shapefile de recorte
  3. use gdal_translatepara cortar até as extensões da forma
  4. use gdalwarpcom -cutlineparâmetro

Os passos 2 e 3 são para otimização, você pode conviver com apenas gdalwarp -cutline ....

Consulte Cortando rasters com GDAL usando polígonos do Linfinity para solução baseada em linux, todos agrupados em um script. Outro exemplo de linha de corte pode ser visto no tutorial de Michael Corey, criando colinas para o Mapnik .

Matt Wilson
fonte
Matt, você deve se lembrar que trac.osgeo.org/gdal/ticket/1599 parece que a linha de recorte cumpre isso
Mike T
10

Parece que esse assunto está sempre voltando. Eu mesmo não sabia que o GDAL> 1.8 é tão avançado que já lhe fornece um tratamento justo da linha de comando para executar essa tarefa.

O comentário de Mike Toews é bastante útil, mas você pode simplesmente fazer, por exemplo:

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

Você pode agrupar esse comando dentro de um script python com o excelente módulo de subprocesso .

Uma coisa realmente problemática para mim é que eu precisava fornecer uma solução mínima para esse problema, o que significa o mais simples possível e não exige muitas dependências externas. O uso da Python Imaging Library, como no tutorial de Joel Lawhead, é puro, mas eu vim com a seguinte solução: usando matrizes mascaradas Numpy.
Não sei se é melhor, mas era o que eu sabia do que (há 3 anos ...).
Originalmente, criei uma área de dados válida dentro da varredura original (por exemplo, a extensão da varredura de saída onde era a mesma), mas gostei da idéia de tornar a varredura também menor (por exemplo, -crop_to_cutline), por isso adotei world2Pixelde Joel Lawhead. Aqui está minha própria solução:

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)

para obter uma descrição completa dos class MaskRastermétodos e seus, consulte o github do meu projeto .

Usando esse código, você ainda precisará usar o GDAL. No entanto, o plano é usar no futuro o Python puro, sempre que possível, porque o público-alvo do meu software tem dificuldades com muitas dependências (eu uso o Debian para desenvolver o software e os clientes usam o Windows 7 ...).

Oz123
fonte
Eu gosto do exemplo da linha de comando que você deu, mas você pode explicar o que o argumento -crop_to_cutline faz? Não sei ao certo qual é o seu objetivo, pois o shapefile de recorte é especificado por -cutline.
23412 hendra
11
a opção -cutline recorta a varredura na caixa delimitadora interna da camada de polígono. Por exemplo, se for menor nas extensões, o raster de saída também será menor. Sem isso, o raster de saída terá o mesmo tamanho do original, mas com NULL em todos os pontos fora da área de seu interesse.
Oz123