Correct way to define projection of Raster from NumPy Array?

Discussion created by amilt on Jul 4, 2013
Latest reply on Jul 7, 2013 by SWalbridge-esristaff
I am attempting to convert a NumPy array to a Raster using arcpy.NumPyArrayToRaster and am having trouble with referencing the Raster to properly place it in space.

My basic steps are this:

  1. Take a raster and project it into a new coordinate system using ProjectRaster_management

  2. Convert the entire Raster to a numpy array using arcpy.RasterToNumPyArray

  3. Take a section of the NumPy array based on a new lower-left corner and number of cells in the new coordinate system

  4. Convert the numpy array back to a raster using arcpy.NumPyArrayToRaster

  5. Define the projection of the Raster and finish

The problem arises in step 4 when I convert back to a raster. The command I use is

arcpy.NumPyArrayToRaster(N.astype(pixelType), lowerLeft, cellSize, cellSize, float(noDataValue))

where lowerLeft is the lower left corner in the new coordinate system. When I do this, the resulting Raster has a lower-left corner that is different from the one given to the command (one in command is "512309.548826205 2085425.58116366" and one in output is "849785.857717518 4616397.45350893").

Why does this happen?

I have noticed that the help for NumPyArrayToRaster says the lower-left corner should be in "map units," but I have yet to figure out whether or not this makes a difference. If it does, I do not see why I would get the wrong output since the projection of the raster and the lower-left corner are in the same coordinate system.

Any ideas how to fix this?

EDIT: This issue may be related to another problem I am having, which is that even though I define the projection of the raster to be the same as the the coordinate system I am working in, when I go to the raster, it will define the coordinate system to be in the default for ArcMap (I think). This does not happen when running from Python alone. This still happens when I set arcpy.env.outputCoordinateSystem to be the same as the system I am working in.

Perhaps ArcGIS is forcing the raster into this other coordinate system, which makes the lower-left corner wrong in the original system. When I then DefineProjection for the raster, it must convert back to the original which is incorrect...?

EDIT 2: I just checked and, yes indeed the raster is being forced into this other coordinate system. I see that the system the raster is being projected into is the output coordinate system set within the ArcMap mxd. So, my new question is

Why does setting arcpy.env.outputCoordinateSystem not change the coordinate system used to project the Raster upon conversion?

--------------------------- TEMPORARY SOLUTION --------------------------
EDIT 3: I've found a temporary, terrible solution.

The basic process is to first figure out what coordinate system the converted numpy array will be projected into, to create a dummy reference point so it is properly referenced in that system, and then to project to the actual output system. This, to me, is a bad solution, especially because it forces me to write the output raster to disk, not to mention the several extra steps. This problem needs to be addressed

Here's some code to go before step 4 above:

# get the coordinate system the raster will be projected into (wtf) dummy = arcpy.NumPyArrayToRaster(N.astype(pixelType)) defaultSystem = arcpy.Describe(dummy).spatialReference  # project the lower-left point into this system dummyPointGeometry = arcpy.PointGeometry(lowerLeft, coordSys) dummyPointGeometryProjected = dummyPointGeometry.projectAs(defaultSystem) dummyPoint = dummyPointGeometryProjected.lastPoint  # get the cell size conversion defaultMPU = defaultSystem.metersPerUnit outputMPU = coordSys.metersPerUnit dummyCellSize = ( defaultMPU / float(outputMPU) ) * cellSize  # convert the NumPy array to a raster arcpy.env.outputCoordinateSystem = defaultSystem # this is important. R = arcpy.NumPyArrayToRaster(N.astype(pixelType), dummyPoint, dummyCellSize, dummyCellSize, float(noDataValue))  # project the raster result = arcpy.ProjectRaster_management(R, outraster, coordSys)