I am trying to run ZonalStatisticsAsTable using python.
I have a feature class in a FGDB and a raster within an SDE database (SQL Server).
If I load these into layers within arcmap, open the python console (within arcmap) and execute:
arcpy.CheckOutExtension("Spatial")
arcpy.sa.ZonalStatisticsAsTable(r"Site Boundary\Site", "Site", "sde_jbacfm.JBAGISEDITOR.River09Q1000D", "TempTable")
I get a table of results.
If I run this code, from a python console outside of arcmap:
import arcpy
arcpy.CheckOutExtension("Spatial")
# In these filepaths I have purposefully replaced the full path with '...' for clarity. This is just for this example
arcpy.env.workspace = r"C:\...\Output.gdb" # This is an empty geodatabase initiated with a WGS84 projection
input_feat = r"C:\...\SitePolygon.gdb\Site Boundary\Site"
raster = r"C:\...\connection.sde\sde_jbacfm.JBAGISEDITOR.River09Q1000D"
output = "TempTable"
arcpy.sa.ZonalStatisticsAsTable(input_feat, "Site", raster, output)
This works but I get an entirely different result for the statistics.
I am trying to work out what is going on? The only thing I can think of is that the projection of each input (and the output database) is different, but I had understood that arcpy would handle different projections automatically?
In any case, I would still expect to get the same result using python through arcmap as I did not reproject anything there either.
first this line should be put in raw format arcpy.env.workspace = "C:\...\Output.gdb"
I wouldn't assume anything about...nor would I rely upon...projections on the fly. You can easily fix that by establishing them explicitly
The line is in raw format; just an error in my editing for this forum.
Projections on the fly or not...it doesnt explain why it would behave differently within arcmap than outside?
The two should be the same (right?), so should both give the same result, right or wrong.
When I see a difference in tool output inside and outside ArcMap (especially raster tools) my first place to check is the geoprocessing environment (extent, cellsize, output coordinate system, etc). If the environment settings are different in the ArcMap document than the environment defaults that are used when you run a python script standalone, you are likely to see different results.
> The only thing I can think of is that the projection of each input (and the output database) is different, but I had understood that arcpy would handle different projections automatically?
Yes, but its automatic choice will depend on your geoprocessing settings and the input datasets. This varies by tool. In this situation, you probably want to explicitly set the spatial reference. Zonal Statistics analysis is area-weighted so you should specify an equal-area (projected) coordinate system rather than just letting arcpy geoprocessing guess based in your inputs.
The easiest way to do do this in Python is to use the EPSG code:
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(102039) # Albers equal area
The env information is exactly the same, apart from the output workspace - however, these are both just file geodatabases with the same parameters.
env is mostly the same because the parameters are not set (i.e. extent, cellsize, output coordinate system are all blank in both arcmap and stand-alone).
I would've thought that in this case though, arcpy would always choose the same parameters for any settings it requires which are not user-provided, as the inputs are exactly the same whether it is arcmap or stand-alone.
You are right:
Setting
arcpy.env.outputCoordinateSystem
ensures the same answer whether from within arcmap or stand alone, however it is a different answer than that gained when run from within arcmap without
arcpy.env.outputCoordinateSystem.
We have decided that the latter case is the right answer, so I'm still trying to figure out what has gone wrong.
Sounds like I had your answer. The ArcMap defaults clearly aren't the same as what you expected.
When I do raster processing I try to project the data (Project Raster, Resample) instead of project on the fly, because then I have full control of how the resampling works.
Explicitly controlling the raster processing environment (cell size, snap raster, etc) makes sure you are comparing apples with apples.
Yes I have it working by projecting the input feature to the same projection as the raster to extract stats from.
It does seem a bit non-intuitive that arcpy would opaquely alter settings from within arcmap and not without.
> It does seem a bit non-intuitive that arcpy would opaquely alter settings from within arcmap and not without.
When running from ArcMap, there are settings the user can set in the geoprocessing environment dialog, and in modelbuilder -- and also if you call your tool from a python script, it may have environment settings before your tool is called. To support this, the tool behavior must be affected by the GP environment set in the application. So if you care (like in this case, with zonal stats) you need to explicitly set the env in your tool. Your script documentation should inform the user which geoprocessing settings your tool will honor.
Seems like I've not fully resolved this issue as it appears that trying to project two of the same feature, but with different coordinate systems, to the same coordinate system gives slightly different results.
Importing two features which are the same but have different SRS's to arcmap correctly displays them (i.e. on top of each other).
Reprojecting both to a 3rd coordinate system outside of arcmap and then importing these new features now shows that they are not in the same position.
I'm wondering what other environment settings are controlling this.