래스터 내에서 주어진 픽셀 값의 최소 경계 범위를 찾으십니까? 전역 범위로 설정됩니다. 이

특정 값을 가진 래스터의 최소 경계 범위를 찾는 방법이 있는지 궁금합니다. 전역 이미지에서 래스터를 자르고 범위는 NoData 영역이 많은 전역 범위로 설정됩니다. 이 래스터에서 NoData 영역을 제거하고 특정 값의 픽셀을 포함하는 영역 만 유지하고 싶습니다. 어떻게해야합니까?

다음은 내 예입니다. value = 1 (파란색 영역)을 추출하고 추가 처리를 위해 전 세계가 아닌 파란색 영역의 범위를 사용하고 싶습니다.

샘플 이미지



답변

문제를 올바르게 이해했다면 null이 아닌 값의 최소 경계 상자를 알고 싶은 것처럼 들립니다. 래스터를 다각형으로 변환하고 관심있는 다각형을 선택한 다음 다시 래스터로 변환 할 수 있습니다. 그런 다음 최소 경계 상자를 제공하는 속성 값을 볼 수 있습니다.


답변

트릭은 값이있는 데이터의 한계를 계산하는 것입니다. 아마도 가장 빠르고 가장 자연스럽고 가장 일반적인 방법은 영역 요약을 사용하는 것입니다. 영역에 모든 비 NoData 셀을 사용하면 X 및 Y 좌표를 포함하는 격자의 최소 및 최대 영역이 전체 범위를 제공합니다.

ESRI는 이러한 계산이 수행되는 방식을 계속 변경합니다. 예를 들어, 좌표 그리드의 내장 표현식은 ArcGIS 8로 삭제되었으며 반환되지 않은 것으로 보입니다. 재미를 위해 여기에 무엇이든지 상관없이 빠르고 간단한 계산이 있습니다.

  1. 다음과 같이 그리드를 자체와 동일시 하여 단일 영역으로 변환

    “My grid”== “내 그리드”

  2. 값이 1 인 상수 그리드를 누적하여 열 인덱스 그리드를 만듭니다. (인덱스에 0으로 시작합니다.) 원하는 경우이 값에 셀 크기를 곱하고 원점의 x 좌표를 추가하여 x 좌표 그리드 ” X “(아래 참조).

  3. 마찬가지로, 값이 64 인 상수 그리드를 누적 하여 행 인덱스 그리드 ( 및 y 좌표 그리드 “Y”)를 만듭니다.

  4. 단계 (1)의 구역 격자를 사용하여 “X”및 “Y”의 구역 최소 및 최대계산하십시오 . 이제 원하는 범위를 갖습니다.

(두 개의 구역 통계표에 표시된 범위는이 그림에서 직사각형으로 표시됩니다. 그리드 “I”는 단계 (1)에서 얻은 구역 그리드입니다.)

더 나아가려면 출력 테이블에서이 4 개의 숫자를 추출하고이를 사용하여 분석 범위를 제한해야합니다. 제한된 분석 범위가있는 원본 그리드를 복사하면 작업이 완료됩니다.


답변

다음 은 Python 도구 상자 (.pyt) 인 ArcGIS 10.1+ 용 @whubers 메서드 버전입니다 .

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)

링크가 조금씩 404로 이동하면 스크립트 가 Github의 코드 숨김 상태 입니다. 이 폴더는 일부 재구성에 적합합니다.


답변

모든 분석 성능을 위해 ArcGIS에는 GIMP 와 같은 전통적인 데스크탑 이미지 편집기에서 찾을 수있는 기본 래스터 조작이 없습니다 . 출력 범위 환경 설정 을 수동으로 무시하지 않는 한 출력 래스터에 대해 입력 래스터와 동일한 분석 범위를 사용하려고 합니다. 이것이 바로 설정하려는 것이 아니라 찾고자하는 것이므로 ArcGIS 방식으로 작업을 수행하고 있습니다.

최선의 노력에도 불구하고 프로그래밍에 의지하지 않고 원하는 계산의 부분 집합을 얻을 수있는 방법을 찾을 수 없었습니다 (래스터-벡터 변환없이 계산 낭비).

대신 김프를 사용하여 색상 별 선택 도구를 사용하여 파란색 영역을 선택한 다음 선택 영역을 반전시키고, 나머지 픽셀을 제거하려면 삭제를 누르고, 선택 영역을 다시 반전시키고, 이미지를 선택 영역으로 자르고 마지막으로 다시 내보냈습니다. PNG. 김프는이를 1 비트 깊이 이미지로 저장했습니다. 결과는 다음과 같습니다.

물론 샘플 이미지에 공간 참조 구성 요소가없고 김프가 공간적으로 인식하지 못하기 때문에 출력 이미지는 샘플 입력만큼 유용합니다. 공간 분석에 사용 하려면 지리 참조 가 필요합니다 .


답변

다음은 SAGA GIS를 사용할 수있는 가능성입니다. http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

SAGA GIS에는 작업을 수행하는 “Crop to Data”모듈 (Grid Tools 모듈 라이브러리)이 있습니다.

그러나이를 위해서는 GDAL 가져 오기 모듈로 Geotif를 가져와 SAGA에서 처리 한 다음 GDAL 내보내기 모듈을 사용하여 Geotif로 다시 내 보내야합니다.

단지는 ArcGIS GP 도구를 사용하여 또 다른 가능성은 사용하여 래스터에서 TIN 구축하는 것입니다 TIN에 래스터를 사용하여 그 경계를 계산 TIN 도메인을 , 그리고 클립 경계 (또는 사용하여 봉투하여 래스터를 다각형으로 기능 봉투 ).