开源GIS-geos实现空间快速连接
其中关于空间计算的只有简单的判断
标题
而要真正实现空间连接，是需要将两个shp文件进行里外循环，例如，以其中一个shp的要素个数为外循环，以另一个shp要素数量为内循环，逐一判断是否存在空间上的关联(包含、被包含。。。)
如此做循环肯定会影响程序执行的速度。因此可以使用Rtree作为索引，加快循环的检索。
话不多说，下面上代码，以“包含”的空间关系做空间连接为例
#-----------------打开OSM数据---------------------------
osm_ds = ogr.Open(osmfile)
osm_lyr=osm_ds.GetLayer(0)
feat=osm_lyr.GetFeature(0)
keys_osm=feat.keys()
spatial=osm_lyr.GetSpatialRef()
count=osm_lyr.GetFeatureCount()
#---打开BUffer-----------
buffer_s=ogr.Open(bufferfile)
buff_lyr=buffer_s.GetLayer(0)
buff_feature=buff_lyr.GetFeature(0)
keys_buff=buff_feature.keys()
#--创建索引--
shp_index=index.Index()
#--shp几何数据导入Rtree--
for i,feat in enumerate(buff_lyr):
# geom=feat.GeometryDef()
geom = feat.GetGeometryRef()
bbox=geom.GetEnvelope()#--left，right，buttom，top
# print("Type of bbox",type(bbox),"：",bbox)
shp_index.insert(i,(bbox[0],bbox[2],bbox[1],bbox[3]))#---idx.insert(0, (left, bottom, right, top))
#--空间连接
# print()
gdal.SetConfigOption( "GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption( "SHAPE_ENCODING", "UTF-8")
Spatial_driver = ogr.GetDriverByName("ESRI Shapefile")
spatial_ds = Spatial_driver.CreateDataSource(matchedfile)
match_lyr = spatial_ds.CreateLayer(‘matched‘,spatial) # 创建缓冲区
#创建字段
# print("osm_keys:",keys_osm)
osm_def = osm_lyr.GetLayerDefn()
# for i,each in enumerate(keys_osm):
# fielddef=osm_def.GetFieldDefn(i)
# if fielddef==None:
# break
# # print(fielddef)
# # print(each)
# fieldDefn = ogr.FieldDefn(each,fielddef.GetType())
# match_lyr.CreateField(fieldDefn)
# print("Buffer_keys:", keys_buff)
keys_match=keys_buff.copy()
keys_match.extend(keys_osm)
# print("matchedfile Fields:", keys_match)
buff_def = buff_lyr.GetLayerDefn()
count_field=0
for i,key in enumerate(keys_match):
field=None
if key in keys_buff:
field=buff_def.GetFieldDefn(i)
count_field+=1
else:
field = osm_def.GetFieldDefn(i-count_field)
if field is None:
break
fieldDef= ogr.FieldDefn(key,ogr.OFTString)
fieldDef.SetWidth(254)
match_lyr.CreateField(fieldDef)
#---添加几何属性和字段属性----------
match_feat = ogr.Feature(match_lyr.GetLayerDefn())
# print("--创建空间连接...")
for i,osm_feat in enumerate(osm_lyr):
# print(i+1,"/%d"%count,end=‘,‘)
geom = osm_feat.geometry().Clone()
match_feat.SetGeometry(geom)
# print("type of geom",type(geom))
#--buffer属性
feat_within=[]
#--在Rtree中寻找相交的部分要素
bbox=geom.GetEnvelope()
feat_id = list(shp_index.intersection((bbox[0],bbox[2],bbox[1],bbox[3])))
# for buffer_feat in buff_lyr: 大连人流手术费用
# if geom.Within(buffer_feat.geometry().Clone()):
# feat_within.append(buffer_feat)
#---被包含的要素--
for idx in feat_id:
feat_temp=buff_lyr.GetFeature(idx)
if geom.Within(feat_temp.geometry().Clone()):
feat_within.append(feat_temp)
if len(feat_within)==1:
for key in keys_match:
if key in keys_osm:
match_feat.SetField(key, str(osm_feat.GetField(key)))
# print(key, "(one):", osm_feat.GetField(key), end=‘,‘)
if key in keys_buff:
# print(key, "(one):", feat_within[0].GetField(key), end=‘,‘)
match_feat.SetField(key, (feat_within[0].GetField(key)))
elif len(feat_within)>1:
data = getMostvalue(feat_within)
for key in keys_match:
if key in keys_osm:
match_feat.SetField(key, str(osm_feat.GetField(key)))
# print(key, "(one of):", osm_feat.GetField(key), end=‘,‘)
if key in keys_buff:
# print(key,"(one of):",feat_within[data].GetField(key),end=‘,‘)
match_feat.SetField(key,str(feat_within[data].GetField(key)))
else:
# --osm属性
for key in keys_match:
if key in keys_osm:
match_feat.SetField(key, str(osm_feat.GetField(key)))
# print(key, "(None):", osm_feat.GetField(key), end=‘,‘)
else:
match_feat.SetField(key,"None")
print(key, "(None):","None", end=‘,‘)
match_lyr.CreateFeature(match_feat)
# print()
feat_within.clear()
#--进度条---
# time.sleep(0.5)
# sys.stdout.write("创建空间连接... %2f%% \r" % (100*(i+1)/count))
# sys.stdout.flush()
# print()
spatial_ds.Destroy()
以上是个人在做的一个项目的程序，简单的实现了快速的空间连接，关于Rtree的教程，可以百度搜索Rtree，使用的库包为geos和rtree。