ITPub博客

首页 > 应用开发 > Python > 开源GIS-geos实现空间快速连接

开源GIS-geos实现空间快速连接

原创 Python 作者:ckxllf 时间:2021-03-29 17:01:28 0 删除 编辑

  开源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。

来自 “ ITPUB博客 ” ,链接:http://blog.itpub.net/69945560/viewspace-2765530/,如需转载,请注明出处,否则将追究法律责任。

请登录后发表评论 登录
全部评论

注册时间:2019-08-16

  • 博文量
    198
  • 访问量
    176994