Keeping spatial reference using arcpy.RasterToNumPyArray?

Question

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 Luke's second method

My tentative method My tentative method

Answer 1

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)

• yosukesabai: Luke, Thank you for the answer. It appears to me that I need to do something hybrid of your two methods? Both method sets extent, spref etc. correct, but fail to array data correctly... First method makes all cell to be nodata. Second method works better, but they dont seem to handle nodata cells in myArray correctly (sorry that i didnt tell my cell has some nodata and i want to keep those intact). Some trial and error made me think that I have to take second approach but arcpy.env.outCoordinateSystem makes output looks reasonable. not much understanding of how so though.
• user2856: You didn't ask about NoData, you asked about spatial reference and cell size. The arcpy.NumPyArrayToRaster method takes a value_to_nodata parameter.
• yosukesabai: Luke, Thank you for edit. I tried your method of providing 5th argument (value_to_nodata), but i still yielded the figure at the top (nodata cell filled with 0 or 255, and nodata_value is not set for the output raster). the only workaround i found was to set env.outputCoordinateSystem before NumPyArrayToRaster, instead of using DefineProjection_management afterward. It doens't make sense why this works, but i just go with my solution. thank you for all the help.

Answer 2

I had some problems getting ArcGIS to handle NoData values correctly with the examples shown here. I extended the example from reomtesensing.io blog (which is more or less similar to the solutions shown here) to better handle NoData.

Apparently ArcGIS (10.1) likes the value -3.40282347e+38 as NoData. So I convert back and forth between numpy NaN and -3.40282347e+38. The code is here:

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)