I am using ArcGIS 10.1, and want to create a new raster based on two preexisting rasters. The RasterToNumPyArray has a good example which I want to adapt.
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")
Problem is that it strips the spatial reference and also cell size. I figured it has to do arcpy.env, but how do I set them based on input raster? I cannot figure it out.
Taking Luke's answer, this is my tentative solution.
Both of Luke's solution set spatial reference, extent and cell size correctly. But the first method did not carry data in the array correctly and output raster is filled with nodata everywhere. His second method works mostly, but where i have big region of nodata, it fills with blocky zeros and 255s. This may have to do with how i handled nodata cells, and i am not quite sure how i was doing it (should be another Q though). I included images of what i am talking about.
#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)
Image showing results. I both case nodata cells are shown yellow.
Luke's second method
My tentative method
Answer
Check out the Describe method.
Something like the following should work.
#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)
OR
#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)
EDIT: The The arcpy.NumPyArrayToRaster method takes a value_to_nodata parameter. Use it like so:
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)
No comments:
Post a Comment