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

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.

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.

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.

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

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

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.

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:

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.

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