Correct way to define projection of Raster from NumPy Array?

[AustinMilt]

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 Raster.save 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)

[ShaunWalbridge]

Austin,

That sounds like quite an annoying bug, I imagine it took some time to track down the hacky solution you ended up with. I'd like to reproduce the bug so we can get an issue on this, as from my reading of the help, NumPyArrayToRaster should respect arcpy.env.outputCoordinateSystem. Could you just share a couple of details: What are the versions of OS and ArcGIS that you're using to produce this issue? Do I need to use a certain projection for the map, or anything that differs from the data will trigger this?

cheers,
Shaun

[AustinMilt]

Austin,

That sounds like quite an annoying bug, I imagine it took some time to track down the hacky solution you ended up with. I'd like to reproduce the bug so we can get an issue on this, as from my reading of the help, NumPyArrayToRaster should respect arcpy.env.outputCoordinateSystem. Could you just share a couple of details:

What are the versions of OS and ArcGIS that you're using to produce this issue?

I am using Windows 7 Enterprise SP1 and ArcGIS 10.1 SP1

Do I need to use a certain projection for the map, or anything that differs from the data will trigger this?

Good question. The map projection that produced this error is NAD_1983_UTM_Zone_17N and I was using a custom projection (.prj attached) inside the Python script. I just tried it with the same map projection but with WGS_1984_UTM_Zone_17N as the script coordinate system and got the same results.

[AustinMilt]

Just to add a tip for those who may be experiencing the same problem: If you're testing a Python Toolbox from within ArcMap, you need to start a new instance of the tool in order to overwrite any previous Map spatial reference settings that you have since changed.

[AustinMilt]

Austin,

That sounds like quite an annoying bug, I imagine it took some time to track down the hacky solution you ended up with. I'd like to reproduce the bug so we can get an issue on this, as from my reading of the help, NumPyArrayToRaster should respect arcpy.env.outputCoordinateSystem. Could you just share a couple of details: What are the versions of OS and ArcGIS that you're using to produce this issue? Do I need to use a certain projection for the map, or anything that differs from the data will trigger this?

cheers,
Shaun

I might also add that I seem to be unable to overwrite the Map environment's output coordinate system with arcpy.env.outputCoordinateSystem when saving/outputing a raster using any of the Spatial Analyst in-memory tools.

For instance, one other workaround solution might be to do something like this:

arcpy.env.outputCoordinateSystem = coordSys
mapOpen = True
try: 

    mapDoc = arcpy.mapping.MapDocument('CURRENT')
    mapFrame = mapDoc.activeDataFrame
    mapSR = mapFrame.spatialReference
    mapFrame.spatialReference = coordSys
    
# no map document open
except: mapOpen = False


# # # Convert the Array to Raster # # #

# convert the NumPyArray to a raster
R = arcpy.NumPyArrayToRaster(N.astype(pixelType), lowerLeft, cellSize, cellSize, float(noDataValue))

# save the output raster
R.save(outraster)

# reset the map coordinate system
if mapOpen: mapFrame.spatialReference = mapSR

But I find that the output raster still has the coordinate system set in the current map document's output coordinate system

[ShaunWalbridge]

Just to add a tip for those who may be experiencing the same problem: If you're testing a Python Toolbox from within ArcMap, you need to start a new instance of the tool in order to overwrite any previous Map spatial reference settings that you have since changed.

One way around this is to force the relevant code to reload each time the toolbox is opened. That's what I generally do when developing a toolbox, then remove them for a release.