Нахождение минимальной степени ограничения заданного значения пикселя в растре?

Мне интересно, есть ли способ найти минимальный ограничивающий экстент для растра с определенным значением. Я вырезал растр из глобального изображения, и экстент установлен как глобальный экстент с большой областью NoData. Я хотел бы удалить область NoData из этого растра и сохранить только область, содержащую пиксели определенного значения. Как я могу это сделать?

Вот мой пример: я хотел бы извлечь значение = 1 (синяя область) и использовать экстент синей области, а не весь мир для дальнейшей обработки.

Если я правильно понял вопрос, это звучит так, как будто вы хотите знать минимальную ограничивающую рамку значений, которые не равны нулю. Возможно, вы можете преобразовать растр в полигоны, выбрать интересующие вас полигоны и затем преобразовать их обратно в растр. Затем вы можете посмотреть значения свойств, которые должны дать вам минимальный ограничивающий прямоугольник.

Хитрость заключается в том, чтобы вычислить пределы данных, которые имеют значения. Пожалуй, самый быстрый, самый естественный и самый общий способ их получения - это зональные сводки: при использовании всех ячеек, не относящихся к NoData, для зоны минимальные и максимальные значения сетки, содержащие координаты X и Y, обеспечат полную протяженность.

ESRI постоянно меняет способы выполнения этих расчетов; например, встроенные выражения для координатных сеток были удалены с помощью ArcGIS 8 и, похоже, не вернулись. Просто для удовольствия, вот набор быстрых, простых вычислений, которые сделают работу, несмотря ни на что.

1. Преобразуйте сетку в одну зону , приравняв ее к себе, как в

   "Моя сетка" == "Моя сетка"

2. Создайте сетку индекса столбца, накапливая поток постоянной сеткой со значением 1. (Индексы начнутся с 0.) При желании умножьте это на размер ячейки и добавьте x-координату начала координат, чтобы получить сетку x-координат " Х "(показано ниже).

3. Аналогично, создайте сетку индекса строки ( а затем сетку с координатами Y "Y") путем накопления потока постоянной сеткой со значением 64.

4. Используйте сетку зоны из шага (1), чтобы вычислить минимальные и максимальные значения зон «X» и «Y» : теперь у вас есть желаемый экстент.

(Экстент, как показано в двух таблицах зональной статистики, изображен прямоугольным контуром на этом рисунке. Сетка «I» - это сетка зоны, полученная на шаге (1).)

Чтобы пойти дальше, вам нужно будет извлечь эти четыре числа из их выходных таблиц и использовать их для ограничения степени анализа. Копирование исходной сетки с ограниченным экстентом анализа завершает задачу.

Вот версия метода @whubers для ArcGIS 10.1+ в виде набора инструментов 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')

Я разработал решение на основе gdal и numpy. Он разбивает растровую матрицу на строки и столбцы и удаляет любую пустую строку / столбец. В этой реализации «пустой» - это что-то меньше 1, и учитываются только одноканальные растры.

(Когда я пишу, я понимаю, что этот подход к сканирующей линии подходит только для изображений с «воротниками» нодаты. Если ваши данные представляют собой острова в морях нулей, пространство между островами также будет уменьшено, сведя все вместе и полностью испортив географическую привязку .)

Бизнес-части (нуждаются в уточнении, не будут работать как есть):

    #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

В полном сценарии:

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)

Скрипт находится в моем тайнике с кодом на Github, если ссылка идет 404-й разом; эти папки созрели для некоторой реорганизации.

Несмотря на всю свою аналитическую мощь, в ArcGIS отсутствуют базовые растровые манипуляции, которые можно найти в традиционных настольных редакторах изображений, таких как GIMP . Ожидается, что вы хотите использовать тот же экстент анализа для выходного растра, что и входной, если вы не переопределите параметр среды Выходной экстент вручную . Так как это именно то, что вы хотите найти, а не установить, ArcGIS делает все возможное.

Несмотря на все мои усилия и не прибегая к программированию, я не мог найти способа получить размер желаемого подмножества изображения (без преобразования растра в вектор, которое в вычислительном отношении расточительно).

Вместо этого я использовал GIMP, чтобы выделить синюю область, используя инструмент выбора по цвету, а затем инвертировал выделение, ударил «Удалить», чтобы удалить остальные пиксели, снова инвертировал выделение, обрезал изображение до выделенного и, наконец, экспортировал его обратно в PNG. GIMP сохранил его как изображение глубиной 1 бит. Результат ниже:

Конечно, поскольку в вашем образце не хватает компонента пространственной привязки, а GIMP не знает пространственно, выходное изображение примерно так же полезно, как и ввод данных в качестве образца. Вам нужно будет привязать его, чтобы он был полезен при пространственном анализе.

Вот одна из возможностей использования SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

В SAGA GIS есть модуль «Crop to Data» (в библиотеке модулей Grid Tools), который выполняет задачу.

Но для этого потребуется импортировать Geotif с помощью модуля импорта GDAL, обработать его в SAGA и, наконец, снова экспортировать его как Geotif с помощью модуля экспорта GDAL.

Другая возможность, использующая только инструменты ArcGIS GP, состоит в том, чтобы создать TIN из вашего растра, используя растр в TIN , вычислить его границу, используя домен TIN , и обрезать свой растр по границе (или его конверт, используя Feature Envelope to Polygon ).