¿Extraer el valor del raster para puntos específicos (con coordenadas X e Y) y escribirlo en una columna de la tabla de atributos?

Question

Tengo puntos específicos en un archivo shp y una lista de rasters en el mismo sistema de coordenadas y quiero extraer los valores de los rásteres para cada punto y colocar estos valores en la columna de la Tabla de Atributos del archivo shp de puntos.

Los rastreos son archivos separados.

¿Es posible pasar de vector a trama?

Archivo Shp pf puntos:

infc = 'samplepoints1.shp'

código sólo un comienzo:

import arcpy

infile = r"C:\Users\anna_pavlenko\Desktop\SILVER\biscuit_dnbr.tif"
#create raster object
ftc = arcpy.sa.Raster(infile)
#rows and columns
ftc.height 
ftc.width 
#cell size
cellSize = ftc.meanCellHeight
print 'The cell  size is ', cellSize, 'm'

rasterList = r"C:\Users\biscuit_dnbr.tif"
rstArray = arcpy.RasterToNumPyArray(rasterList)

rows, cols = rstArray.shape
print "number of rows: " ,  rows
print "number of columns: ", cols

for rowNum in xrange(rows):
    for colNum in xrange(cols):
        value = rstArray.item(rowNum, colNum)

Si uso

inRasterList = [["biscuit_dnbr.tif", "foothills_dnbr.tif", "dNBR1"] ,     ["biscuit_rdnbr.tif", "foothills_rdnbr.tif", "RdNBR1"],["silver_dnbr.tif", "pony_complex_dnbr.tif", "dNBR2"], ["silver_rdnbr.tif", "pony_complex_rdnbr.tif", "RdNBR2"]] 
#  ["foothills_dnbr.tif","dNBR1"], "foothills_rdnbr.tif", "RdNBR1"],
#  ["pony_complex_dnbr.tif", "dNBR2"], "pony_complex_rdnbr.tif", "RdNBR2"]]

ExtractMultiValuesToPoints(inPointFeatures, inRasterList, "NONE")

Esta función crea valores para 1 lugar (polígono) y extrae el valor del punto dentro de este polígono solamente y los coloca en columnas: "dNBR1", "dNBR2", "RdNBR1", "RdNBR2".
Para otro conjunto de rásters crea 0 -zeros para los primeros puntos y añade nuevas columnas : "dNBR1_1", "dNBR2_1", "RdNBR1_1", "RdNBR2_" cada vez.

Answer 1

Si no tiene la extensión del analista espacial o quiere utilizar Numpy, puede convertir las coordenadas de los puntos x e y en la fila y la columna que estén cerca de ese punto. Las herramientas de Extracción de Valores tienen muchas más opciones, y producen resultados ligeramente diferentes a los del código de abajo. No estoy defendiendo un método sobre otro. Tampoco sé si esto sería más rápido.

import numpy as np

Obtenga información sobre su archivo raster:

rast = arcpy.Raster(r"somePath")
desc = arcpy.Describe(rast)

Obtenga la esquina superior izquierda de la trama y el sistema de referencia espacial (puede ignorar la referencia espacial si todos sus puntos y tramas son iguales).

ulx = desc.Extent.XMin
uly = desc.Extent.YMax
#Spatialrefernce
sr = desc.spatialReference

Convertir el raster en un array numpy

rstArray  = arcpy.RasterToNumPyArray(rast)

Acceda a su Feature Class utilizando un cursor de actualización para añadir el valor. Utilizo la herramienta projectAs para asegurarme de que el punto está en el mismo sistema de coordenadas que el raster. Calcule la fila y la columna en función de la distancia del punto a la esquina superior izquierda. Esto puede depender mucho del redondeo. A continuación, pasar esto en su matriz numpy.

with arcpy.da.UpdateCursor("NHDPoint",["SHAPE@","value"]) as uc:
    for row in uc:
        pnt = row[0].projectAs(sr)
        #assuming every point falls to the left and below uperleft corner
        deltaX = pnt.centroid.X - ulx
        deltaY = lly- pnt.centroid.Y
        arow = int(deltaY/rast.meanCellHeight)
        acol = int(deltaX/rast.meanCellWidth)
        row[1] = rstArray[arow,acol]
        uc.updateRow(row)

Este código asume que los puntos caerán a la derecha de la esquina superior izquierda, y por debajo de ella. También acaba de utilizar el redondeo por defecto al convertir el flotador a un entero.

Otra opción, que podría estar más cerca de lo que las herramientas de Extracto es utilizar una distancia al centro de la celda de sus puntos para determinar qué valor utilizar.

Answer 2

Además de la herramienta señalado por HDunn puede utilizar Extraer los multivalores a los puntos si tiene muchos datos ráster con el mismo shapefile de puntos. Aquí hay otro fragmento de la ayuda de ArcGIS.

import arcpy
from arcpy.sa import *
from arcpy import env 
env.workspace = "c:/sapyexamples/data"
ExtractMultiValuesToPoints("observers.shp", [["elevation", "ELEV"], 
                           ["costraster", "COST"], ["flowdir", "DIR"]], "NONE")

Answer 3

Si ya estás usando ArcPy, puedes usar el Extract Values to Points herramienta (suponiendo que tenga la extensión Spatial Analyst)

Este es el fragmento de código del sitio de ESRI

ExtractValuesToPoints("rec_sites.shp", "elevation",
                  "C:/sapyexamples/output/outValPnts","INTERPOLATE",
                  "VALUE_ONLY")

Iterar sobre sus capas rasterizadas como lo hizo y extraer los valores a su capa de puntos