Before responding to your points... let me clarify...
My code works... but I don't like the answer.
When I spatially identify the research sites that conform to a set of rules, ~25% of the selections don't meet the criteria. I know this because the selected table items contain the actual grid values and 25% of the field property values clearly don't match the criteria.
Because the file contains the actual grid values associated with each research site, I can apply the rules directly to the shapefile table. And when I do that... ~25% of the selections (though true) aren't selected when I do the spatial analysis.
Varying grid cell size is the primary source of the spatial selection error. ArcGIS rounds cells to the coarsest size / grid orientation before calculating the math (in this case... inference space rules).
I can resolve the problem by changing the cell size / grid orientation of all the grids to match that of the highest resolution cell. The ~1km cell would become 16 ~250m cells and the ~500m cell would become 4 ~250m cells. That's option a.
I can convert the highest resolution grid to points (file geodatabase), extract the 7 grid values associated with each point, and apply the inference space rules to the extracted values in the file geodatabase. The selected points represent the entire inference space. Then I'd need to use the selected points to identify the associated research site points. One approach: convert the selected points back into a ~250m grid, convert the grid to polygons and intersect the polygon with the research sites. That's option b.
Option c: Convert the highest resolution grid to points (file geodatabase), extract the 7 grid values associated with each point and convert those points into an integer grid. (The grid format might not be the best raster format for this approach due to the number of cells and table values.)
Crop production is only relevant in that I'll need an areal representation of all locations that meet the inference space criteria (whether or not there are research sites at those cells) in order to characterize the production potential of that area. "All locations" means all of Sub-Saharan Africa.
The table also contains crop x nutrient modeling coefficients developed from research site data. Some coefficients are based on a considerable amount of research. Others are not. The modeling coefficients will allow a researcher to locate potential undeveloped production areas.
As for your comments...
>The extract to point tool extracts a grid value for a point and creates a new point file. But the Extract Multi Values to Points tool modifies (i.e. adds the extracted point values) to the original point file.
However... I did run into a data corruption issue. I usually create point shape files from tab-delimited text files exported from Excel. The coefficient values are all floating point. But... if the first row of coefficients were integers ArcGIS turned all the coefficient data into integers. Now I'm creating point files from Excel files (no extra rows or columns) and I haven't run into any problems.
>TableSelect worked for me. Here's the code for selecting extracted dem values from the research site point file (Data9-28-15.shp) using the "dem" value associated with a user-selected point:
# ======= 6 of 7 =========
demmax = 5818
print "dem500 resolution ~500m; Range: -155 - 5818; Selected dem500 value: " + str(dem)
# Decrease amount to add (demAdd) if dem + 300 is undefined
if demmax - dem < 300:
demAdd = demmax - dem
##print "Will add " + str(demAdd) + " instead of 300"
else:
demAdd = 300
##print "Will add 300"
demGT = dem - 300
demLT = dem + demAdd
demField = arcpy.AddFieldDelimiters("Data9-28-15.shp","DEM500")
# Create rule-based where clause for TableSelect
if dem < 700:
whereClause = "{0} <= {1}".format(demField,1000)
print whereClause
else:
whereClause = "{0} > {1} AND {0} < {2}".format(demField,demGT,demLT)
print whereClause
# Select data and output to table
arcpy.TableSelect_analysis("D:/New/.../Out/sandout", "D:/New/.../Out/demout", whereClause)
My arcpy/Python background consists of the 2-3 hour ESRI online short courses (Python Basics & Python Scripting for ArcGIS), a book and a whole lot of internet searching. I don't know anything about python/numpy
I think I've addressed your other points.
