Trouver l'étendue limite minimale d'une valeur de pixel donnée dans le raster?

Je me demande s'il existe un moyen de trouver l'étendue de délimitation minimale pour un raster avec une valeur particulière. J'ai découpé un raster à partir d'une image globale et l'étendue est définie comme l'étendue globale avec beaucoup de zone NoData. Je voudrais supprimer la zone NoData de ce raster et ne conserver que la zone contenant les pixels de la valeur particulière. Comment puis-je faire ceci?

Voici mon exemple: je voudrais extraire la valeur = 1 (zone bleue) et utiliser l'étendue de la zone bleue plutôt que le monde entier pour un traitement ultérieur.

Si j'ai bien compris la question, il semble que vous souhaitiez connaître la zone de délimitation minimale des valeurs qui ne sont pas nulles. Vous pouvez peut-être convertir le raster en polygones, sélectionner les polygones qui vous intéressent, puis les reconvertir en raster. Vous pouvez ensuite regarder les valeurs des propriétés qui devraient vous donner la boîte englobante minium.

L'astuce consiste à calculer les limites des données qui ont des valeurs. Le résumé peut être le moyen le plus rapide, le plus naturel et le plus général d'obtenir ces résumés: en utilisant toutes les cellules non NoData pour la zone, les min et max zonaux des grilles contenant les coordonnées X et Y fourniront l'étendue complète.

ESRI ne cesse de changer la façon dont ces calculs peuvent être effectués; par exemple, les expressions intégrées pour les grilles de coordonnées ont été supprimées avec ArcGIS 8 et semblent ne pas être revenues. Juste pour le plaisir, voici un ensemble de calculs simples et rapides qui feront le travail quoi qu'il arrive.

1. Convertissez la grille en une seule zone en l'assimilant à elle-même, comme dans

   "Ma grille" == "Ma grille"

2. Créez une grille d'index de colonne en accumulant en flux une grille constante de valeur 1. (Les index commenceront par 0.) Si vous le souhaitez, multipliez-la par la taille de la cellule et ajoutez la coordonnée x de l'origine pour obtenir une grille de coordonnées x " X "(illustré ci-dessous).

3. De même, créez une grille d'index de ligne ( puis une grille de coordonnées y "Y") en accumulant en flux une grille constante de valeur 64.

4. Utilisez la grille de zone de l'étape (1) pour calculer le min et le max zonal de "X" et "Y" : vous avez maintenant l'étendue souhaitée.

(L'étendue, comme indiqué dans les deux tableaux de statistiques zonales, est représentée par un contour rectangulaire sur cette figure. La grille "I" est la grille de zone obtenue à l'étape (1).)

Pour aller plus loin, vous devrez extraire ces quatre nombres de leurs tableaux de sortie et les utiliser pour limiter l'étendue de l'analyse. La copie de la grille d'origine, avec l'étendue d'analyse restreinte en place, termine la tâche.

Voici une version de la méthode @whubers pour ArcGIS 10.1+ en tant que boîte à outils 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')

J'ai conçu une solution basée sur gdal et numpy. Il divise la matrice raster en lignes et colonnes et supprime toute ligne / colonne vide. Dans cette implémentation, "vide" est inférieur à 1 et seuls les rasters à bande unique sont pris en compte.

(Je me rends compte au moment où j'écris que cette approche de la ligne de balayage ne convient qu'aux images avec des "colliers" nodata. Si vos données sont des îles dans des mers nulles, l'espace entre les îles sera également supprimé, écrasant tout ensemble et gâchant totalement le géoréférencement .)

Les parties commerciales (doivent être développées, ne fonctionneront pas telles quelles):

    #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

Dans un script complet:

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)

Le script est dans ma cachette de code sur Github, si le lien va 404 chasser un peu; ces dossiers sont mûrs pour une réorganisation.

Malgré toute sa puissance analytique, ArcGIS manque de manipulations raster de base que vous pouvez trouver avec les éditeurs d'images de bureau traditionnels comme GIMP . Il s'attend à ce que vous souhaitiez utiliser la même étendue d'analyse pour votre raster en sortie que votre raster en entrée, sauf si vous remplacez manuellement le paramètre d'environnement Output Extent . Étant donné que c'est exactement ce que vous cherchez à trouver, et non à définir, la façon de faire d'ArcGIS se met en travers.

Malgré mes meilleurs efforts, et sans recourir à la programmation, je n'ai trouvé aucun moyen d'obtenir l'étendue de votre sous-ensemble souhaité de l'image (sans conversion raster en vecteur qui est un gaspillage de calcul).

Au lieu de cela, j'ai utilisé GIMP pour sélectionner la zone bleue à l'aide de l'outil de sélection par couleur, puis inversé la sélection, appuyez sur Supprimer pour supprimer le reste des pixels, inversé à nouveau la sélection, recadré l'image à la sélection et enfin exporté vers PNG. GIMP l'a enregistré en tant qu'image de profondeur 1 bit. Le résultat est ci-dessous:

Bien sûr, parce que votre échantillon d'image ne comportait pas de composant de référence spatiale et que GIMP n'est pas conscient de l'espace, l'image de sortie est à peu près aussi utile que votre entrée d'échantillon. Vous devrez le géoréférencer pour qu'il soit utile dans une analyse spatiale.

Voici une possibilité en utilisant SAGA GIS: http://permalink.gmane.org/gmane.comp.gis.gdal.devel/33021

Dans SAGA GIS, il existe un module "Crop to Data" (dans la bibliothèque de modules Grid Tools), qui effectue la tâche.

Mais cela vous obligerait à importer votre Geotif avec le module d'importation GDAL, à le traiter dans SAGA, et enfin à l'exporter à nouveau en tant que Geotif avec le module d'exportation GDAL.

Une autre possibilité en utilisant uniquement les outils ArcGIS GP serait de créer un TIN à partir de votre raster à l'aide de Raster vers TIN , de calculer sa frontière à l'aide du domaine TIN et de découper votre raster par la frontière (ou son enveloppe à l'aide de l'enveloppe d' entité au polygone ).