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
Mi metodo tentativo
fuente
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í:
fuente