#! /usr/bin/env pythonimportsysimportfionafromshapely.geometryimportshape,mappingfromshapely.wktimportdumps,loadsfile=sys.argv[1]target=sys.argv[2]expression=sys.argv[3]d=fiona.open(file)outputSchema={"geometry":d.schema["geometry"],"properties":{}}fields=[]geomField=Noneforiinrange(4,len(sys.argv)):fieldExpression=sys.argv[i]# field parsingasIndex=fieldExpression.find(" as ")fieldNameAndType=fieldExpression[:asIndex].strip()fieldEvalExpression=fieldExpression[asIndex+4:]colonIndex=fieldNameAndType.find(":")ifcolonIndex!=-1:fieldName=fieldNameAndType[:colonIndex]fieldType=fieldNameAndType[colonIndex+1:]computed=Trueelse:fieldName=fieldNameAndTypefieldType=d.schema["properties"][fieldEvalExpression]computed=Falsefield={"name":fieldName,"type":fieldType,"expression":fieldEvalExpression,"computed":computed}iffield["name"]=="geometry":geomField=fieldelse:fields.append(field)# create field in new featureoutputSchema["properties"][field["name"]]=field["type"]ifgeomFieldisnotNone:outputSchema["geometry"]=geomField["type"]iflen(fields)==0:outputSchema["properties"]=d.schema["properties"]output=fiona.open(target,"w",driver="ESRI Shapefile",crs=d.crs,schema=outputSchema)forfeatureind:ifeval(expression):newFeature={"geometry":feature["geometry"],"properties":{}}ifgeomFieldisnotNone:newFeature["geometry"]=eval(geomField["expression"])# If there are no field ops include alliflen(fields)==0:newFeature["properties"]=feature["properties"]else:forfieldinfields:# field evaluationvalue=Noneiffield["computed"]:# Expressionvalue=eval(field["expression"])else:# Just field referencevalue=feature["properties"][field["expression"]]# create field in new featurenewFeature["properties"][field["name"]]=valueoutput.write(newFeature)d.close()output.close()

#! /usr/bin/env pythonimportsysimportcollectionsimportfionafromshapely.geometryimportshape,mapping,MultiPoint,MultiLineString,MultiPolygonfromshapely.opsimportcascaded_unionfile=sys.argv[1]target=sys.argv[2]geometryType=sys.argv[3]d=fiona.open(file)outputSchema={"geometry":geometryType,"properties":{}}groupField=NonegroupFieldUsed=Noneiflen(sys.argv)>4:groupField=sys.argv[4]groupFieldUsed=TrueoutputSchema["properties"][groupField]=d.schema["properties"][groupField]else:groupField="id"groupFieldUsed=FalseoutputSchema["properties"]["id"]="int"output=fiona.open(target,"w",driver="ESRI Shapefile",crs=d.crs,schema=outputSchema)classes={}counter=0total=len(d)forfeatureind:print"\rgroup:\t",50*counter/total,counter=counter+1ifgroupFieldUsed:value=feature["properties"][groupField]else:value=0ifvalueinclasses:class_=classes[value]else:class_=[]classes[value]=class_class_.append(feature)counter=0total=len(classes)forvalueinclasses:print"\rgroup:\t",50+50*counter/total,counter=counter+1class_=classes[value]classGeometries=[shape(feature["geometry"])forfeatureinclass_]unionResult=cascaded_union(classGeometries)# hack because cascaded_union may not give a collectionifnotisinstance(unionResult,collections.Iterable):ifgeometryType=="MultiPoint":unionResult=MultiPoint([unionResult])elifgeometryType=="MultiLineString":unionResult=MultiLineString([unionResult])elifgeometryType=="MultiPolygon":unionResult=MultiPolygon([unionResult])feature={"geometry":mapping(unionResult),"properties":{groupField:value}}output.write(feature)d.close()output.close()

defgetField(fieldExpression,schema):asIndex=fieldExpression.find(" as ")fieldNameAndType=fieldExpression[:asIndex].strip()fieldEvalExpression=fieldExpression[asIndex+4:]colonIndex=fieldNameAndType.find(":")ifcolonIndex!=-1:fieldName=fieldNameAndType[:colonIndex]fieldType=fieldNameAndType[colonIndex+1:]computed=Trueelse:fieldName=fieldNameAndTypefieldType=schema["properties"][fieldEvalExpression]computed=Falsereturn{"name":fieldName,"type":fieldType,"expression":fieldEvalExpression,"computed":computed}defgetFields(args,schema):fields=[]forfieldExpressioninargs:fields.append(getField(fieldExpression,schema))returnfieldsdefgetGeometryField(fields):returnnext((fieldforfieldinfieldsiffield["name"]=="geometry"),None)defgetAlphanumericFields(fields):return[fieldforfieldinfieldsiffield["name"]!="geometry"]

Y con el siguiente script (shape_join.py):

#! /usr/bin/env pythonimportschema_parserimportsysimporttimeimportfionafromshapely.geometryimportshape,mappingfromshapely.opsimportcascaded_unionfromrtreeimportindexclassSequentialScan:defpreLoop(self):passdeffeaturesFor(self,outerFeature,inner):returninnerclassSpatialIndexScan:innerMemory=[]idx=index.Index()defpreLoop(self,inner):# Load inner in memory for random accessforinnerFeatureininner:self.innerMemory.append(innerFeature)bounds=shape(innerFeature["geometry"]).boundsself.idx.insert(len(self.innerMemory)-1,bounds)deffeaturesFor(self,outerFeature,inner):ret=[]# Query the indexqueryResult=self.idx.intersection(shape(outerFeature["geometry"]).bounds)forinnerFeatureIndexinqueryResult:ret.append(self.innerMemory[innerFeatureIndex])returnretouterPath=sys.argv[1]innerPath=sys.argv[2]target=sys.argv[3]scanType=sys.argv[4]joinCondition=sys.argv[5]start=time.time()outer=fiona.open(outerPath)inner=fiona.open(innerPath)ifscanType=="sequential":innerScan=SequentialScan()elifscanType=="spatial-index":innerScan=SpatialIndexScan()innerScan.preLoop(inner)combinedSchema=dict(outer.schema.items()+inner.schema.items())fields=schema_parser.getFields(sys.argv[6:],combinedSchema)iflen(fields)==0:print"field expressions missing"sys.exit(-1)else:outputSchema={"properties":{}}geomField=schema_parser.getGeometryField(fields)ifgeomFieldisNone:print"geometry field expression missing"sys.exit(-1)else:outputSchema["geometry"]=geomField["type"]alphanumericFields=schema_parser.getAlphanumericFields(fields)forfieldinalphanumericFields:outputSchema["properties"][field["name"]]=field["type"]output=fiona.open(target,"w",driver="ESRI Shapefile",crs=outer.crs,schema=outputSchema)counter=0total=len(outer)forouterFeatureinouter:print"\rjoin:\t\t",100*counter/total,counter=counter+1scannedFeatures=innerScan.featuresFor(outerFeature,inner)forinnerFeatureinscannedFeatures:ifeval(joinCondition):newFeature={"geometry":eval(geomField["expression"]),"properties":{}}forfieldinalphanumericFields:# field evaluationvalue=eval(field["expression"])# create field in new featurenewFeature["properties"][field["name"]]=valueoutput.write(newFeature)output.close()inner.close()outer.close()end=time.time()printend-start,"seconds"