Encontrando a extensão limite mínima de um determinado valor de pixel na varredura?

9

Gostaria de saber se existe uma maneira de encontrar a extensão limite mínima para uma varredura com um valor específico. Recortei uma varredura de uma imagem global e a extensão é definida como extensão global com muita área NoData. Gostaria de remover a área NoData desta varredura e reter apenas a área que contém os pixels do valor específico. Como posso fazer isso?

Aqui está o meu exemplo: eu gostaria de extrair valor = 1 (área azul) e usar a extensão da área azul em vez de o mundo inteiro para processamento adicional.

Imagem de exemplo

Visto
fonte
Você poderia postar uma amostra?
Aaron
"Gostaria de excluir as linhas e colunas nulas desta varredura." O que exatamente isso significa? Não entendo qual é o resultado final desejado.
precisa saber é o seguinte
Por "extensão limite mínima", você está procurando uma extensão retangular ou uma pegada poligonal representando a área da imagem com dados?
blah238
11
@ Tomome, o OP está procurando descobrir a extensão, não precisa criá-lo manualmente.
blah238
11
Se literalmente alguma coisa é justa, então algum software possui comandos internos para fazer isso; consulte reference.wolfram.com/mathematica/ref/ImageCrop.html, por exemplo.
whuber

Respostas:

6

Se eu entendi a pergunta corretamente, parece que você deseja saber a caixa delimitadora mínima dos valores que não são nulos.Talvez você possa converter a varredura em polígonos, selecione os polígonos nos quais está interessado e depois converta-os novamente em uma varredura. Você pode ver os valores das propriedades que devem fornecer a caixa delimitadora mínima.

dango
fonte
11
No total, essa é provavelmente a melhor abordagem, dados os limites do fluxo de trabalho de processamento raster do ArcGIS.
blah238
Eu fiz exatamente isso. Existe alguma maneira automática? Eu acho que o algoritmo de varredura para polígono tem um passo para extrair a caixa delimitadora mínima da varredura.
Visto
você está atrás de uma solução python?
Dango #
8

O truque é calcular os limites dos dados que possuem valores. Talvez a maneira mais rápida, mais natural e mais geral de obtê-las seja com resumos zonais: usando todas as células não NoData para a zona, o zonal min e max das grades contendo as coordenadas X e Y fornecerão toda a extensão.

A ESRI continua mudando as maneiras pelas quais esses cálculos podem ser feitos; por exemplo, expressões integradas para as grades de coordenadas foram descartadas com o ArcGIS 8 e parecem não ter retornado. Apenas por diversão, aqui está um conjunto de cálculos rápidos e simples que farão o trabalho, não importa o quê.

  1. Converta a grade em uma única zona , equiparando-a a si mesma, como em

    "Minha grade" == "Minha grade"

  2. Crie uma grade de índice de coluna acumulando uma grade constante com o valor 1. (Os índices começarão com 0.) Se desejar, multiplique pelo tamanho da célula e adicione a coordenada x da origem para obter uma grade com coordenadas x " X "(mostrado abaixo).

  3. Da mesma forma, crie uma grade de índice de linha ( e, em seguida, uma grade de coordenadas y "Y") acumulando uma grade constante com o valor 64.

  4. Use a grade de zona da etapa (1) para calcular o zonal mínimo e máximo de "X" e "Y" : agora você tem a extensão desejada.

Imagem final

(A extensão, conforme mostrado nas duas tabelas de estatísticas zonais, é representada por um contorno retangular nesta figura. Grade "I" é a grade de zona obtida na etapa (1).)

Para ir além, você precisará extrair esses quatro números de suas tabelas de saída e usá-los para limitar a extensão da análise. A cópia da grade original, com a extensão de análise restrita, conclui a tarefa.

whuber
fonte
8

Aqui está uma versão do método @whubers para ArcGIS 10.1+ como uma caixa de ferramentas python (.pyt).

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Raster Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [ClipNoData]


class ClipNoData(object):
    def __init__(self):
        """Clip raster extent to the data that have values"""
        self.label = "Clip NoData"
        self.description = "Clip raster extent to the data that have values. "
        self.description += "Method by Bill Huber - https://gis.stackexchange.com/a/55150/2856"

        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        params = []

        # First parameter
        params+=[arcpy.Parameter(
            displayName="Input Raster",
            name="in_raster",
            datatype='GPRasterLayer',
            parameterType="Required",
            direction="Input")
        ]

        # Second parameter
        params+=[arcpy.Parameter(
            displayName="Output Raster",
            name="out_raster",
            datatype="DERasterDataset",
            parameterType="Required",
            direction="Output")
        ]

        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return arcpy.CheckExtension('spatial')==u'Available'

    def execute(self, parameters, messages):
        """See https://gis.stackexchange.com/a/55150/2856
           ##Code comments paraphrased from @whubers GIS StackExchange answer
        """
        try:
            #Setup
            arcpy.CheckOutExtension('spatial')
            from arcpy.sa import *
            in_raster = parameters[0].valueAsText
            out_raster = parameters[1].valueAsText

            dsc=arcpy.Describe(in_raster)
            xmin=dsc.extent.XMin
            ymin=dsc.extent.YMin
            mx=dsc.meanCellWidth
            my=dsc.meanCellHeight
            arcpy.env.extent=in_raster
            arcpy.env.cellSize=in_raster
            arcpy.AddMessage(out_raster)

            ## 1. Convert the grid into a single zone by equating it with itself
            arcpy.AddMessage(r'1. Convert the grid into a single zone by equating it with itself...')
            zones = Raster(in_raster) == Raster(in_raster)

            ## 2. Create a column index grid by flow-accumulating a constant grid
            ##    with value 1. (The indexes will start with 0.) Multiply this by
            ##    the cellsize and add the x-coordinate of the origin to obtain
            ##    an x-coordinate grid.
            arcpy.AddMessage(r'Create a constant grid...')
            const = CreateConstantRaster(1)

            arcpy.AddMessage(r'2. Create an x-coordinate grid...')
            xmap = (FlowAccumulation(const)) * mx + xmin

            ## 3. Similarly, create a y-coordinate grid by flow-accumulating a
            ##    constant grid with value 64.
            arcpy.AddMessage(r'3. Create a y-coordinate grid...')
            ymap = (FlowAccumulation(const * 64)) * my + ymin

            ## 4. Use the zone grid from step (1) to compute the zonal min and
            ##    max of "X" and "Y"
            arcpy.AddMessage(r'4. Use the zone grid from step (1) to compute the zonal min and max of "X" and "Y"...')

            xminmax=ZonalStatisticsAsTable(zones, "value", xmap,r"IN_MEMORY\xrange", "NODATA", "MIN_MAX")
            xmin,xmax=arcpy.da.SearchCursor(r"IN_MEMORY\xrange", ["MIN","MAX"]).next()

            yminmax=ZonalStatisticsAsTable(zones, "value", ymap,r"IN_MEMORY\yrange", "NODATA", "MIN_MAX")
            ymin,ymax=arcpy.da.SearchCursor(r"IN_MEMORY\yrange", ["MIN","MAX"]).next()

            arcpy.Delete_management(r"IN_MEMORY\xrange")
            arcpy.Delete_management(r"IN_MEMORY\yrange")

            # Write out the reduced raster
            arcpy.env.extent = arcpy.Extent(xmin,ymin,xmax+mx,ymax+my)
            output = Raster(in_raster) * 1
            output.save(out_raster)

        except:raise
        finally:arcpy.CheckInExtension('spatial')
user2856
fonte
Muito bom Luke. Independente, executável, usa in_memory e bem comentado para inicializar. Eu tive que desativar o processamento em segundo plano ( Geoprocessamento> opções> ... ) para fazê-lo funcionar.
Matt Wilkie
11
Atualizei o script e configurei canRunInBackground = False. Observarei que vale a pena alterar os ambientes de espaço de trabalho / scratchworkspace para uma pasta local (não FGDB), como encontrei quando os deixei como padrão (ou seja, <perfil de usuário de rede> \ Default.gdb), o script levou 9min !!! para executar em uma grade de 250x250 células. Mudar para um FGDB locais levou 9sec e para uma pasta local 1-2sec ...
user2856
bom argumento sobre pastas locais e obrigado pela rápida correção em segundo plano (muito melhor do que escrever instruções para todos os que forneço). Pode valer a pena vomitar isso no bitbucket / github / gcode / etc.
Matt Wilkie
+1 Obrigado por esta contribuição, Luke! Agradeço por preencher a lacuna (bastante grande) deixada na minha resposta.
whuber
4

Eu inventei uma solução baseada em gdal e numpy. Ele divide a matriz de varredura em linhas e colunas e elimina qualquer linha / coluna vazia. Nesta implementação, "vazio" é menor que 1, e apenas os rasters de banda única são responsáveis.

(Percebo enquanto escrevo que essa abordagem da linha de varredura é adequada apenas para imagens com "colares" nodata. Se seus dados forem ilhas em mares de nulos, o espaço entre as ilhas também será reduzido, juntando tudo e estragando totalmente o georreferenciamento. .)

As partes do negócio (as necessidades precisam ser aprimoradas, não funcionarão como estão):

    #read raster into a numpy array
    data = np.array(gdal.Open(src_raster).ReadAsArray())
    #scan for data
    non_empty_columns = np.where(data.max(axis=0)>0)[0]
    non_empty_rows = np.where(data.max(axis=1)>0)[0]
        # assumes data is any value greater than zero
    crop_box = (min(non_empty_rows), max(non_empty_rows),
        min(non_empty_columns), max(non_empty_columns))

    # retrieve source geo reference info
    georef = raster.GetGeoTransform()
    xmin, ymax = georef[0], georef[3]
    xcell, ycell = georef[1], georef[5]

    # Calculate cropped geo referencing
    new_xmin = xmin + (xcell * crop_box[0]) + xcell
    new_ymax = ymax + (ycell * crop_box[2]) - ycell
    cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

    # crop
    new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

    # write to disk
    band = out_raster.GetRasterBand(1)
    band.WriteArray(new_data)
    band.FlushCache()
    out_raster = None

Em um script completo:

import os
import sys
import numpy as np
from osgeo import gdal

if len(sys.argv) < 2:
    print '\n{} [infile] [outfile]'.format(os.path.basename(sys.argv[0]))
    sys.exit(1)

src_raster = sys.argv[1]
out_raster = sys.argv[2]

def main(src_raster):
    raster = gdal.Open(src_raster)

    # Read georeferencing, oriented from top-left
    # ref:GDAL Tutorial, Getting Dataset Information
    georef = raster.GetGeoTransform()
    print '\nSource raster (geo units):'
    xmin, ymax = georef[0], georef[3]
    xcell, ycell = georef[1], georef[5]
    cols, rows = raster.RasterYSize, raster.RasterXSize
    print '  Origin (top left): {:10}, {:10}'.format(xmin, ymax)
    print '  Pixel size (x,-y): {:10}, {:10}'.format(xcell, ycell)
    print '  Columns, rows    : {:10}, {:10}'.format(cols, rows)

    # Transfer to numpy and scan for data
    # oriented from bottom-left
    data = np.array(raster.ReadAsArray())
    non_empty_columns = np.where(data.max(axis=0)>0)[0]
    non_empty_rows = np.where(data.max(axis=1)>0)[0]
    crop_box = (min(non_empty_rows), max(non_empty_rows),
        min(non_empty_columns), max(non_empty_columns))

    # Calculate cropped geo referencing
    new_xmin = xmin + (xcell * crop_box[0]) + xcell
    new_ymax = ymax + (ycell * crop_box[2]) - ycell
    cropped_transform = new_xmin, xcell, 0.0, new_ymax, 0.0, ycell

    # crop
    new_data = data[crop_box[0]:crop_box[1]+1, crop_box[2]:crop_box[3]+1]

    new_rows, new_cols = new_data.shape # note: inverted relative to geo units
    #print cropped_transform

    print '\nCrop box (pixel units):', crop_box
    print '  Stripped columns : {:10}'.format(cols - new_cols)
    print '  Stripped rows    : {:10}'.format(rows - new_rows)

    print '\nCropped raster (geo units):'
    print '  Origin (top left): {:10}, {:10}'.format(new_xmin, new_ymax)
    print '  Columns, rows    : {:10}, {:10}'.format(new_cols, new_rows)

    raster = None
    return new_data, cropped_transform


def write_raster(template, array, transform, filename):
    '''Create a new raster from an array.

        template = raster dataset to copy projection info from
        array = numpy array of a raster
        transform = geo referencing (x,y origin and pixel dimensions)
        filename = path to output image (will be overwritten)
    '''
    template = gdal.Open(template)
    driver = template.GetDriver()
    rows,cols = array.shape
    out_raster = driver.Create(filename, cols, rows, gdal.GDT_Byte)
    out_raster.SetGeoTransform(transform)
    out_raster.SetProjection(template.GetProjection())
    band = out_raster.GetRasterBand(1)
    band.WriteArray(array)
    band.FlushCache()
    out_raster = None
    template = None

if __name__ == '__main__':
    cropped_raster, cropped_transform = main(src_raster)
    write_raster(src_raster, cropped_raster, cropped_transform, out_raster)

O script está no meu repositório de códigos no Github, se o link for procurar um pouco; essas pastas estão prontas para alguma reorganização.

Matt Wilson
fonte
11
Isso funcionará para conjuntos de dados realmente grandes. Eu recebo umMemoryError Source raster (geo units): Origin (top left): 2519950.0001220703, 5520150.00012207 Pixel size (x,-y): 100.0, -100.0 Columns, rows : 42000, 43200 Traceback (most recent call last): File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 72, in <module> cropped_raster, cropped_transform = main(src_raster) File "D:/11202067_COACCH/local_checkout/crop_raster.py", line 22, in main data = np.array(raster.ReadAsArray()) MemoryError
user32882
2

Por todo o seu poder analítico, o ArcGIS não possui manipulações básicas de varredura que você pode encontrar com os editores de imagens de desktop tradicionais, como o GIMP . Ele espera que você queira usar a mesma extensão de análise para sua varredura de saída que sua varredura de entrada, a menos que substitua manualmente a configuração do ambiente Extensão de Saída . Como é exatamente isso que você procura, não definido, a maneira de fazer as coisas no ArcGIS está atrapalhando.

Apesar dos meus melhores esforços, e sem recorrer à programação, não consegui encontrar a extensão do subconjunto desejado da imagem (sem conversão de varredura em vetor que é computacionalmente um desperdício).

Em vez disso, usei o GIMP para selecionar a área azul usando a ferramenta Selecionar por cor e inverti a seleção, pressione Excluir para remover o restante dos pixels, inverti a seleção novamente, recortei a imagem para seleção e finalmente a exportei novamente para PNG. O GIMP salvou-o como uma imagem de profundidade de 1 bit. O resultado está abaixo:

Resultado

Obviamente, como sua imagem de amostra não possuía um componente de referência espacial e o GIMP não possui reconhecimento espacial, a imagem de saída é tão útil quanto a entrada de amostra. Você precisará georeferenciar para que seja útil em uma análise espacial.

blah238
fonte
Na verdade, essa operação costumava ser fácil nas versões anteriores do Spatial Analyst: o máximo e o mínimo zonal das duas grades de coordenadas (X e Y), usando o recurso como zona, fornecem exatamente a extensão. (Bem, convém expandi-lo pela metade do tamanho da célula nas quatro direções.) No ArcGIS 10, você precisa ser criativo (ou usar o Python) para criar uma grade de coordenadas. Independentemente disso, tudo pode ser feito inteiramente no SA usando apenas operações de grade e nenhuma intervenção manual.
whuber
@whuber Vi sua solução em outro lugar, mas ainda não tenho certeza de como posso implementar seu método. Você poderia me mostrar mais alguns detalhes disso?
Visto
@Seen Uma rápida pesquisa neste site encontra uma conta do método em gis.stackexchange.com/a/13467 .
whuber
1

Aqui está uma possibilidade usando o SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

No SAGA GIS, existe um módulo "Crop to Data" (na biblioteca do módulo Grid Tools), que executa a tarefa.

Mas isso exigiria que você importasse seu Geotif com o módulo de importação GDAL, processasse-o no SAGA e, finalmente, exportasse como Geotif novamente com o módulo de exportação GDAL.

Outra possibilidade usando apenas as ferramentas do ArcGIS GP seria criar um NIF a partir do seu raster usando Raster para TIN , calcular seu limite usando o domínio TIN e Recortar o seu raster pelo limite (ou seu envelope usando o recurso Envelope para o polígono ).

blah238
fonte