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,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 ) / tmp ) * 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