Mantener referencia espacial usando arcpy.RasterToNumPyArray?

13

Estoy usando ArcGIS 10.1 y quiero crear un nuevo ráster basado en dos rásteres preexistentes. El RasterToNumPyArray tiene un buen ejemplo que quiero adaptarse.

import arcpy
import numpy 
myArray = arcpy.RasterToNumPyArray('C:/data/inRaster')
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum
newRaster = arcpy.NumPyArrayToRaster(myArrayPerc)
newRaster.save("C:/output/fgdb.gdb/PercentRaster")

El problema es que elimina la referencia espacial y también el tamaño de la celda. Pensé que tenía que hacer arcpy.env, pero ¿cómo los configuro en función del ráster de entrada? Yo no puedo averiguarlo.


Tomando la respuesta de Luke, esta es mi solución tentativa.

Ambas soluciones de Luke establecen la referencia espacial, la extensión y el tamaño de celda correctamente. Pero el primer método no transportaba datos en la matriz correctamente y el ráster de salida está lleno de nodos en todas partes. Su segundo método funciona principalmente, pero donde tengo una gran región de nodata, se llena con ceros en bloque y 255. Esto puede tener que ver con cómo manejé las células de nodata, y no estoy muy seguro de cómo lo estaba haciendo (aunque debería ser otra Q). Incluí imágenes de lo que estoy hablando.

#Setting the raster properties directly 
import arcpy 
import numpy 

inRaster0='C:/workspace/test0.tif' 
inRaster1='C:/workspace/test1.tif' 
outRaster='C:/workspace/test2.tif' 

dsc=arcpy.Describe(inRaster0) 
sr=dsc.SpatialReference 
ext=dsc.Extent 
ll=arcpy.Point(ext.XMin,ext.YMin) 

# sorry that i modify calculation from my original Q.  
# This is what I really wanted to do, taking two uint8 rasters, calculate 
# the ratio, express the results as percentage and then save it as uint8 raster.
tmp = [ np.ma.masked_greater(arcpy.RasterToNumPyArray(_), 100) for _ in inRaster0, inRaster1]
tmp = [ np.ma.masked_array(_, dtype=np.float32) for _ in tmp]
tmp = ((tmp[1] ) / tmp[0] ) * 100
tmp = np.ma.array(tmp, dtype=np.uint8)
# i actually am not sure how to properly carry the nodata back to raster...  
# but that's another Q
tmp = np.ma.filled(tmp, 255)

# without this, nodata cell may be filled with zero or 255?
arcpy.env.outCoordinateSystem = sr

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight) 

newRaster.save(outRaster) 

Imagen que muestra resultados. Las dos celdas de nodata del caso se muestran amarillas.

El segundo método de Luke El segundo método de Luke

Mi metodo tentativo Mi metodo tentativo

yosukesabai
fuente

Respuestas:

15

Echa un vistazo al método Describe .

Algo como lo siguiente debería funcionar.

#Using arcpy.env
import arcpy
import numpy

inRaster='C:/workspace/test1.tif'
outRaster='C:/workspace/test2.tif'

dsc=arcpy.Describe(inRaster)
arcpy.env.extent=dsc.Extent
arcpy.env.outputCoordinateSystem=dsc.SpatialReference
arcpy.env.cellSize=dsc.meanCellWidth

myArray = arcpy.RasterToNumPyArray(r)
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc)
newRaster.save(outRaster)

O

#Setting the raster properties directly
import arcpy
import numpy

inRaster='C:/workspace/test1.tif'
outRaster='C:/workspace/test2.tif'

dsc=arcpy.Describe(inRaster)
sr=dsc.SpatialReference
ext=dsc.Extent
ll=arcpy.Point(ext.XMin,ext.YMin)

myArray = arcpy.RasterToNumPyArray(inRaster)
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight)
arcpy.DefineProjection_management(newRaster, sr)

newRaster.save(outRaster)

EDITAR: El método The arcpy.NumPyArrayToRaster toma un parámetro value_to_nodata. Úselo así:

try:
    noDataValue=dsc.noDataValue
    arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight,noDataValue)
except AttributeError: #no data is not defined
    arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight)
usuario2856
fuente
Luke, gracias por la respuesta. ¿Me parece que necesito hacer algo híbrido de sus dos métodos? Ambos métodos establecen la extensión, spref, etc., correctos, pero no pueden agrupar los datos correctamente ... El primer método hace que todas las celdas sean nodata. El segundo método funciona mejor, pero parece que no manejan las células de nodata en myArray correctamente (lo siento, no dije que mi célula tiene algunos nodata y quiero mantenerlos intactos). Algunas pruebas y errores me hicieron pensar que tengo que tomar un segundo enfoque, pero arcpy.env.outCoordinateSystem hace que la salida se vea razonable. Sin embargo, no hay mucha comprensión de cómo.
yosukesabai
1
No preguntaste por NoData, preguntaste por referencia espacial y tamaño de celda. El método arcpy.NumPyArrayToRaster toma un parámetro value_to_nodata.
usuario2856
Luke, gracias por editar. Intenté su método de proporcionar el quinto argumento (value_to_nodata), pero aún así arrojé la figura en la parte superior (la celda de nodata se llena con 0 o 255, y nodata_value no está configurado para el ráster de salida). la única solución que encontré fue establecer env.outputCoordinateSystem antes de NumPyArrayToRaster, en lugar de usar DefineProjection_management después. No tiene sentido por qué esto funciona, pero simplemente voy con mi solución. Gracias por toda la ayuda.
yosukesabai
1

Tuve algunos problemas para que ArcGIS manejara los valores NoData correctamente con los ejemplos que se muestran aquí. Extendí el ejemplo del blog reomtesensing.io (que es más o menos similar a las soluciones que se muestran aquí) para manejar mejor NoData.

Aparentemente, a ArcGIS (10.1) le gusta el valor -3.40282347e + 38 como NoData. Así que convierto de ida y vuelta entre NaN numpy y -3.40282347e + 38. El código está aquí:

import arcpy
import numpy as np

infile  = r'C:\data.gdb\test_in'
outfile = r'C:\data.gdb\test_out'

# read raster with No Data as numpy NaN
in_arr  = arcpy.RasterToNumPyArray(infile,nodata_to_value = np.nan)

# processing
out_arr = in_arr * 10

# convert numpy NaN to -3.40282347e+38
out_arr[np.isnan(out_arr)] = -3.40282347e+38

# information on input raster
spatialref = arcpy.Describe(infile).spatialReference
cellsize1  = arcpy.Describe(infile).meanCellHeight
cellsize2  = arcpy.Describe(infile).meanCellWidth
extent     = arcpy.Describe(infile).Extent
pnt        = arcpy.Point(extent.XMin,extent.YMin)

# save raster
out_ras = arcpy.NumPyArrayToRaster(out_arr,pnt,cellsize1,cellsize2, -3.40282347e+38)
out_ras.save(outfile)
arcpy.DefineProjection_management(outfile, spatialref)
Mathiaskopo
fuente