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

9

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.

Przykładowy obraz

Widziany
źródło
Czy możesz opublikować próbkę?
Aaron
„Chciałbym usunąć puste wiersze i kolumny dla tego rastra”. Co to dokładnie znaczy? Nie rozumiem, jaki jest pożądany efekt końcowy.
blah238
Czy poprzez „minimalny zakres ograniczający” szukasz zasięgu prostokątnego lub wielokąta reprezentującego obszar obrazu z danymi?
blah238
1
@Tomek, PO stara się dowiedzieć , w jakim stopniu, nie trzeba go utworzyć ręcznie.
blah238
1
Jeśli dosłownie wszystko jest uczciwą grą, niektóre programy mają wbudowane polecenia, aby to zrobić; patrz na przykład reference.wolfram.com/mathematica/ref/ImageCrop.html .
whuber

Odpowiedzi:

6

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

Dango
źródło
1
Wszystko wskazuje na to, że jest to prawdopodobnie najlepsze podejście, biorąc pod uwagę przepływ pracy przetwarzania rastrowego ArcGIS.
blah238
Zrobiłem to dokładnie. Czy jest jakiś automatyczny sposób? Myślę, że algorytm raster do wielokąta ma krok w celu wyodrębnienia minimalnej ramki granicznej rastra.
Widziany
szukasz rozwiązania python?
dango
8

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.

Ostateczny obraz

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

Whuber
źródło
8

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')
użytkownik2856
źródło
Bardzo miły Luke. Samodzielny, uruchamialny, używa in_memory i dobrze skomentowany do uruchomienia. Musiałem wyłączyć przetwarzanie w tle ( Geoprocessing> opcje> ... ), aby działało.
matt wilkie
1
Zaktualizowałem skrypt i ustawiłem canRunInBackground = False. Zwrócę uwagę, że warto zmienić środowiska workspace / scratchworkspace na folder lokalny (nie FGDB), co znalazłem, gdy pozostawiłem je jako domyślne (tj. <Profil użytkownika sieciowego> \ Default.gdb) skrypt zajął 9 minut !!! do działania na siatce komórek 250 x 250. Przejście na lokalny FGDB zajęło 9 sekund, a do lokalnego folderu 1-2sek ...
user2856
dobra uwaga na temat folderów lokalnych i dziękuję za szybką poprawkę w tle (znacznie lepiej niż pisanie instrukcji dla wszystkich, którym ją przekazuję). Może warto rzucić to na bitbucket / github / gcode / itp.
matt wilkie
+1 Dzięki za ten wkład, Luke! Doceniam wypełnianie (dość dużej) luki pozostałej w mojej odpowiedzi.
whuber
4

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.

matowe wilkie
źródło
1
Nie zadziała to w przypadku naprawdę dużych zestawów danych. DostajęMemoryError 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

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:

Wynik

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.

blah238
źródło
Rzeczywiście, ta operacja stosuje się łatwo w poprzednich wersjach przestrzennej Analityka: strefowe maksymalna i minimalna z dwóch siatek współrzędnych (X i Y), przy użyciu funkcji jako strefa, otrzymując dokładnie stopniu. (Cóż, możesz chcieć powiększyć go o połowę wielkości komórki we wszystkich czterech kierunkach). W ArcGIS 10 musisz być kreatywny (lub używać Pythona), aby utworzyć siatkę współrzędnych. Niezależnie od tego, wszystko można zrobić całkowicie w ramach SA przy użyciu tylko operacji siatki i bez ręcznej interwencji.
whuber
@ whuber Widziałem twoje rozwiązanie gdzie indziej, ale wciąż nie jestem pewien, jak mogę wdrożyć twoją metodę. Czy możesz mi pokazać więcej szczegółów na ten temat?
Widziany
@Seen Szybkie wyszukiwanie w tej witrynie pozwala znaleźć opis metody pod adresem gis.stackexchange.com/a/13467 .
whuber
1

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

blah238
źródło