Tuesday 22 December 2015

raster - Keeping spatial reference using arcpy.RasterToNumPyArray?


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



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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...