Znalezienie minimalnego zakresu granicznego danej wartości piksela w rastrze?

Zastanawiam się, czy istnieje sposób na znalezienie minimalnego zasięgu granicznego dla rastra o określonej wartości. Wyciągnąłem raster z obrazu globalnego, a zasięg jest ustawiony jako zasięg globalny z dużą ilością obszaru NoData. Chciałbym usunąć obszar NoData z tego rastra i zachować tylko obszar zawierający piksele o określonej wartości. W jaki sposób mogę to zrobić?

Oto mój przykład: Chciałbym wyodrębnić wartość = 1 (niebieski obszar) i użyć zasięgu niebieskiego obszaru zamiast całego świata do dalszego przetwarzania.

JEŻELI poprawnie zrozumiałem pytanie, brzmi to tak, jakbyś chciał znać minimalną ramkę ograniczającą wartości, które nie są zerowe. Być może możesz przekonwertować raster na wielokąty, wybierz interesujące Cię wielokąty, a następnie przekonwertować je z powrotem na raster. Następnie możesz spojrzeć na wartości właściwości, które powinny dać ci minimalną ramkę ograniczającą.

Sztuką jest obliczenie limitów danych, które mają wartości. Być może najszybszym, najbardziej naturalnym i najogólniejszym sposobem ich uzyskania jest zastosowanie streszczeń strefowych: przy użyciu wszystkich komórek spoza strefy NoData dla strefy, minimalne i maksymalne strefy siatki zawierające współrzędne X i Y zapewnią pełny zakres.

ESRI nieustannie zmienia sposoby wykonywania tych obliczeń; na przykład wbudowane wyrażenia dla siatek współrzędnych zostały usunięte z ArcGIS 8 i wydaje się, że nie powróciły. Dla zabawy, oto zestaw szybkich, prostych obliczeń, które wykonają zadanie bez względu na wszystko.

1. Przekształć siatkę w jedną strefę , zrównując ją ze sobą, jak w

   „My grid” == „My grid”

2. Utwórz siatkę indeksu kolumny, gromadząc przepływ stałą siatkę o wartości 1. (Indeksy zaczną się od 0.) W razie potrzeby pomnóż to przez wielkość komórki i dodaj współrzędną x źródła, aby uzyskać siatkę współrzędnych x ” X ”(pokazano poniżej).

3. Podobnie utwórz siatkę indeksu wierszy ( a następnie siatkę współrzędnych Y „Y”), akumulując przepływ stałą siatkę o wartości 64.

4. Użyj siatki stref z kroku (1), aby obliczyć strefową min. I maks. „X” i „Y” : masz teraz pożądany zasięg.

(Zasięg, jak pokazano w dwóch tabelach statystyk strefowych, jest przedstawiony na tym rysunku prostokątnym konturem. Siatka „I” jest siatką strefy uzyskaną w kroku (1).)

Aby przejść dalej, musisz wyodrębnić te cztery liczby z ich tabel wyjściowych i użyć ich do ograniczenia zakresu analizy. Kopiowanie oryginalnej siatki, z ograniczonym zakresem analizy, kończy zadanie.

Oto wersja metody @whubers dla ArcGIS 10.1+ jako przybornika 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')

Opracowałem rozwiązanie oparte na gdal i numpy. Dzieli macierz rastrową na wiersze i kolumny i upuszcza puste wiersze / kolumny. W tej implementacji „pusta” jest mniejsza niż 1 i uwzględniono tylko rastry z pojedynczym pasmem.

(Zdaję sobie sprawę, gdy piszę, że to podejście ze skanowaniem jest odpowiednie tylko dla obrazów z „kołnierzami” nodata. Jeśli twoje dane są wyspami na morzach zerowych, przestrzeń między wyspami również zostanie usunięta, zgniatając wszystko razem i całkowicie psując georeferencje .)

Części biznesowe (wymagają rozwinięcia, nie będą działać tak, jak są):

    #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

W pełnym skrypcie:

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)

Skrypt znajduje się w mojej skrytce kodu na Githubie, jeśli link zawiera 404 polowanie; te foldery są gotowe do reorganizacji.

Mimo całej swojej mocy analitycznej ArcGIS nie ma podstawowych manipulacji rastrowych, które można znaleźć w tradycyjnych edytorach graficznych, takich jak GIMP . Oczekuje, że chcesz użyć tego samego zakresu analizy dla rastra wyjściowego, co raster wejściowy, chyba że ręcznie zastąpisz ustawienie środowiska Output Extent . Ponieważ jest to dokładnie to, czego szukasz, a nie ustawienie, przeszkadza ci sposób ArcGIS.

Pomimo moich najlepszych starań i bez uciekania się do programowania, nie mogłem znaleźć sposobu na uzyskanie pożądanego podzbioru obrazu (bez konwersji rastra na wektor, która jest marnotrawstwem obliczeniowym).

Zamiast tego użyłem GIMP do wybrania niebieskiego obszaru za pomocą narzędzia Wybierz według koloru, a następnie odwróciłem zaznaczenie, uderzyłem w Usuń, aby usunąć resztę pikseli, ponownie odwróciłem zaznaczenie, przyciąłem obraz do zaznaczenia, a na koniec wyeksportowałem go z powrotem do PNG. GIMP zapisał go jako 1-bitową głębię obrazu. Wynik jest poniżej:

Oczywiście, ponieważ w twoim obrazie próbnym brakowało komponentu odniesienia przestrzennego, a GIMP nie jest przestrzennie świadomy, obraz wyjściowy jest tak samo użyteczny jak twój przykładowy plik wejściowy. Będziesz musiał georeferencję , aby mogła być użyta w analizie przestrzennej.

Oto jedna z możliwości korzystania z SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

W SAGA GIS znajduje się moduł „Kadrowanie do danych” (w bibliotece modułów Grid Tools), który wykonuje to zadanie.

Wymagałoby to jednak zaimportowania Geotifu za pomocą modułu importu GDAL, przetworzenia go w SAGA i wreszcie wyeksportowania go jako Geotif za pomocą modułu eksportu GDAL.

Inną możliwością używając tylko narzędzi ArcGIS GP byłoby zbudować TIN korzystając ze swojego rastra Raster Tin , obliczyć jego granica użyciu TIN domeny i Klip swój raster przez granicę (lub koperty przy użyciu Feature Envelope wielokąt ).