arcgis-10.0pythonspatial-join
I am trying to run the Spatial Join function through the Python shell in ArcGIS 10. When you do the point-and-click method of spatial join, you are given a tool window with an option that asks "How do you want the attributes to be summarized?" However, the Python syntax doesn't show a parameter for this option.
How do I control the attribute summarization through Python? Why is there no parameter for this?
Fiona returns Python dictionaries and you can not use poly['properties']['score']) += point['properties']['score']) with a dictionary.
Example of summing attributes using the references given by Mike T:
# read the shapefiles
import fiona
from shapely.geometry import shape
polygons = [pol for pol in fiona.open('poly.shp')]
points = [pt for pt in fiona.open('point.shp')]
# attributes of the polygons
for poly in polygons:
print poly['properties']
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])
OrderedDict([(u'score', 0)])
# attributes of the points
for pt in points:
print i['properties']
OrderedDict([(u'score', 1)])
.... # (same for the 8 points)
Now, we can use two methods, with or without a spatial index:
# iterate through points
for i, pt in enumerate(points):
point = shape(pt['geometry'])
#iterate through polygons
for j, poly in enumerate(polygons):
if point.within(shape(poly['geometry'])):
# sum of attributes values
polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']
# Create the R-tree index and store the features in it (bounding box)
from rtree import index
idx = index.Index()
for pos, poly in enumerate(polygons):
idx.insert(pos, shape(poly['geometry']).bounds)
#iterate through points
for i,pt in enumerate(points):
point = shape(pt['geometry'])
# iterate through spatial index
for j in idx.intersection(point.coords[0]):
if point.within(shape(polygons[j]['geometry'])):
polygons[j]['properties']['score'] = polygons[j]['properties']['score'] + points[i]['properties']['score']
Result with the two solutions:
for poly in polygons:
print poly['properties']
OrderedDict([(u'score', 2)]) # 2 points in the polygon
OrderedDict([(u'score', 1)]) # 1 point in the polygon
OrderedDict([(u'score', 1)]) # 1 point in the polygon
What is the difference ?
After:
schema = fiona.open('poly.shp').schema
with fiona.open ('output.shp', 'w', 'ESRI Shapefile', schema) as output:
for poly in polygons:
output.write(poly)
To go further, look at Using Rtree Spatial Indexing With OGR, Shapely, Fiona
@FelixIP answered this in the comments...
The spatial join is unnecessary when there is a summary statistics function (although the spatial join with 'have their centres' in as the match option should work, it did not for me).
This is useful when there are large polygons with a stat, say, a population estimate. You need to display this data on a common grid (vector). Sometimes two polygons or more intersect with a single grid poly. The following can be used:
Best Answer
The optional
field_mappingparameter is what you're looking for. The page you linked to links to "Mapping input fields to output fields", which will get you started. ThemergeRuleproperty of theFieldMapobjects is the one that controls which type of statistic/summary to calculate.Basically you have to create a single
FieldMappingsobject, which is a collection of individualFieldMapobjects, each with its ownmergeRule.