Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

EDIT:

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

share|improve this question
add comment

1 Answer 1

up vote 7 down vote accepted

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)
share|improve this answer
    
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. –  yosukesabai Oct 22 '12 at 15:31
1  
You didn't ask about NoData, you asked about spatial reference and cell size. The arcpy.NumPyArrayToRaster method takes a value_to_nodata parameter. –  Luke Oct 23 '12 at 2:58
    
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. –  yosukesabai Oct 25 '12 at 16:23
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.