Menemukan batas batas minimum dari nilai piksel yang diberikan dalam raster?

Saya bertanya-tanya apakah ada cara untuk menemukan batas batas minimum untuk raster dengan nilai tertentu. Saya memotong raster dari gambar global dan sejauh ini ditetapkan sebagai tingkat global dengan banyak area NoData. Saya ingin menghapus area NoData dari raster ini dan hanya mempertahankan area yang mengandung piksel dari nilai tertentu. Bagaimana saya bisa melakukan ini?

Berikut adalah contoh saya: Saya ingin mengekstraksi nilai = 1 (Area biru) dan menggunakan luas area biru daripada seluruh dunia untuk diproses lebih lanjut.

JIKA saya telah memahami pertanyaan dengan benar, sepertinya Anda ingin tahu kotak batas minimum dari nilai-nilai yang tidak nol. Mungkin Anda bisa mengonversi raster ke poligon, pilih poligon yang Anda minati dan kemudian mengubahnya kembali menjadi raster. Anda kemudian dapat melihat nilai properti yang akan memberi Anda kotak pembatas minium.

Caranya adalah dengan menghitung batas data yang memiliki nilai. Mungkin cara tercepat, paling alami, dan paling umum untuk memperoleh ini adalah dengan ringkasan zonal: dengan menggunakan semua sel non-NoData untuk zona tersebut, zonal min dan maks kisi-kisi yang berisi koordinat X dan Y akan memberikan cakupan penuh.

ESRI terus mengubah cara penghitungan ini dapat dilakukan; misalnya, ekspresi bawaan untuk kisi-kisi koordinat dijatuhkan dengan ArcGIS 8 dan tampaknya tidak kembali. Hanya untuk bersenang-senang, inilah serangkaian perhitungan cepat dan sederhana yang akan melakukan pekerjaan apa pun yang terjadi.

1. Konversi kisi menjadi satu zona dengan menyamakannya dengan dirinya sendiri, seperti pada

   "Kisi saya" == "Kisi saya"

2. Buat kisi indeks kolom dengan mengumpulkan aliran kisi konstan dengan nilai 1. (Indeks akan mulai dengan 0.) Jika diinginkan, kalikan ini dengan cellsize dan tambahkan koordinat x asal untuk mendapatkan kisi koordinat x " X "(ditampilkan di bawah).

3. Demikian pula, buat kisi indeks baris ( dan kemudian kisi koordinat y "Y") dengan mengakumulasi aliran kisi konstan dengan nilai 64.

4. Gunakan kisi zona dari langkah (1) untuk menghitung zonal min dan maks "X" dan "Y" : Anda sekarang memiliki tingkat yang diinginkan.

(Luasnya, seperti yang ditunjukkan dalam dua tabel statistik zonal, digambarkan oleh garis persegi panjang pada gambar ini. Kotak "I" adalah kotak zona yang diperoleh pada langkah (1).)

Untuk melangkah lebih jauh, Anda perlu mengekstraksi empat angka ini dari tabel output mereka dan menggunakannya untuk membatasi tingkat analisis. Menyalin kisi asli, dengan batas analisis terbatas di tempat, menyelesaikan tugas.

Berikut adalah versi metode @whubers untuk ArcGIS 10.1+ sebagai python toolbox (.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')

Saya telah menemukan solusi berbasis gdal dan numpy. Ini memecah matriks raster menjadi baris dan kolom dan menjatuhkan baris / kolom kosong. Dalam implementasi ini "kosong" adalah kurang dari 1, dan hanya raster band tunggal yang diperhitungkan.

(Saya menyadari ketika saya menulis bahwa pendekatan scanline ini hanya cocok untuk gambar dengan "kerah" nodata. Jika data Anda adalah pulau-pulau di lautan nol, ruang antar pulau akan turun juga, menyelipkan semuanya bersama-sama, dan benar-benar mengacaukan georeferencing .)

Bagian bisnis (perlu disempurnakan, tidak akan berfungsi apa adanya):

    #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

Dalam skrip lengkap:

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)

Skripnya ada di simpanan kode saya di Github, jika tautannya berbunyi 404; folder ini sudah matang untuk reorganisasi.

Untuk semua kekuatan analitisnya, ArcGIS tidak memiliki manipulasi raster dasar yang dapat Anda temukan dengan editor gambar desktop tradisional seperti GIMP . Itu mengharapkan bahwa Anda ingin menggunakan tingkat analisis yang sama untuk raster output Anda sebagai raster input Anda, kecuali Anda secara manual mengabaikan pengaturan lingkungan Extent Output . Karena ini adalah persis apa yang Anda cari, bukan ditetapkan, cara ArcGIS dalam melakukan sesuatu menghalangi.

Terlepas dari upaya terbaik saya, dan tanpa menggunakan pemrograman, saya tidak dapat menemukan cara untuk mendapatkan tingkat subset gambar yang Anda inginkan (tanpa konversi raster-ke-vektor yang secara komputasi boros).

Sebagai gantinya, saya menggunakan GIMP untuk memilih area biru menggunakan alat pilih oleh warna dan kemudian membalikkan seleksi, tekan Delete untuk menghapus sisa piksel, membalikkan seleksi lagi, memotong gambar untuk seleksi, dan akhirnya mengekspornya kembali ke PNG. GIMP menyimpannya sebagai gambar kedalaman 1-bit. Hasilnya di bawah ini:

Tentu saja, karena gambar sampel Anda tidak memiliki komponen referensi spasial, dan GIMP tidak sadar secara spasial, gambar keluaran kira-kira sama berguna dengan input sampel Anda. Anda perlu melakukan georeferensi agar dapat digunakan dalam analisis spasial.

Berikut adalah satu kemungkinan menggunakan SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

Dalam SAGA GIS ada modul "Pangkas ke Data" (di pustaka modul Grid Tools), yang melakukan tugas.

Tetapi ini akan mengharuskan Anda untuk mengimpor Geotif Anda dengan modul impor GDAL, memprosesnya di SAGA, dan akhirnya mengekspornya sebagai Geotif lagi dengan modul ekspor GDAL.

Kemungkinan lain hanya menggunakan alat ArcGIS GP adalah membangun TIN dari raster Anda menggunakan Raster ke TIN , menghitung batasnya menggunakan Domain TIN , dan Klip raster Anda dengan batas (atau amplopnya menggunakan Fitur Envelope to Polygon ).