開源GIS-geos實現空間快速連線

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