¿Manteniendo referencias espaciales 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 adaptar.

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 celda. Me di cuenta de que tiene que hacer arcpy.env, pero ¿cómo los configuro en función del ráster de entrada? No puedo entenderlo.

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

Ambas soluciones de Luke configuran correctamente la referencia espacial, la extensión y el tamaño de la celda. Pero el primer método no llevó datos en la matriz correctamente y el ráster de salida está lleno de nodata en todas partes. Su segundo método funciona principalmente, pero donde tengo una gran región de nodata, se llena con ceros en bloques y 255s. Esto puede tener que ver con la forma en que manejé las células nodatas, y no estoy muy seguro de cómo lo estaba haciendo (aunque debería haber 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. En ambos casos las células nodatas se muestran en amarillo.

El segundo método de Luke

Mimétodotentativo

    
pregunta yosukesabai 22.10.2012 - 05:41

2 respuestas

15

Consulte el 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 arcpy.NumPyArrayToRaster toma un valor_to_nodata parámetro. Úsalo 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)
    
respondido por el Luke 22.10.2012 - 07:01
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 un lado a otro entre números NaN 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)
    
respondido por el Mathiaskopo 01.03.2015 - 19:11

Lea otras preguntas en las etiquetas