from osgeo import ogr
import pandas as pd
defextract_single_ring( point_list_str ):
ring_instance = ogr.Geometry(ogr.wkbLinearRing)
point_list = point_list_str.split(',')
count =0
first_x =None
first_y =Nonefor item in point_list:# 空格分割
str_list = item.split(' ')if(len(str_list)==2):
ring_instance.AddPoint_2D(float(str_list[0]),float(str_list[1]))else:print(str_list[0])if(count ==0):
first_x =float(str_list[0])
first_y =float(str_list[1])
count = count+1
ring_instance.AddPoint_2D(first_x, first_y)return ring_instance
"""
Now I have a list of ring objects. I want to group these rings by area.
In one area, one ring at least intersects with another. please use python.
finally, I want to get the grouped ring with one list that contains lists of rings.
"""defgroup_rings(ogr_ring_list):
grouped_rings =[]while(ogr_ring_list):
current_ring = ogr_ring_list.pop(0)
group =[current_ring]
i =0while i <len(ogr_ring_list):if current_ring.Intersect(ogr_ring_list[i]):
group.append(ogr_ring_list.pop(i))else:
i +=1
grouped_rings.append(group)return grouped_rings
if __name__ =='__main__':# Create a GeoJSON driver
driver = ogr.GetDriverByName("GeoJSON")# Create a new GeoJSON data source
data_source = driver.CreateDataSource("my_geojson_file.geojson")# Create a new layer within the data source with a specific geometry type (e.g., Point)
layer = data_source.CreateLayer("my_layer", geom_type=ogr.wkbUnknown)# Define the layer's fields (e.g., name)
layer.CreateField(ogr.FieldDefn("id", ogr.OFTString))
layer.CreateField(ogr.FieldDefn("ext_path", ogr.OFTString))
layer.CreateField(ogr.FieldDefn("pid", ogr.OFTString))
layer.CreateField(ogr.FieldDefn("deep", ogr.OFTString))
layer.CreateField(ogr.FieldDefn("name", ogr.OFTString))
layer.CreateField(ogr.FieldDefn("geo_wkt", ogr.OFTString))map={'id':0,'pid':1,'deep':2,'name':3,'ext_path':4,'geo':5,'polygon':6}
df = pd.read_csv('ok_geo.csv',encoding='utf-8')#### LOOP ####for index, row in df.iterrows():if(index ==0):continue
geo_wkt_str =None
feature = ogr.Feature(layer.GetLayerDefn())
id_str =str(row.loc['id'])
feature.SetField("id", id_str)
feature.SetField("ext_path",str(row.loc['ext_path']))
feature.SetField("pid",str(row.iloc[map['pid']]))
feature.SetField("deep",str(row.iloc[map['deep']]))
feature.SetField("name",str(row.iloc[map['name']]))
feature.SetField("geo_wkt",f"POINT ({row.iloc[map['geo']]})")
poly_str = row.iloc[map['polygon']]
center_str = row.iloc[map['geo']]# 關注POLY邏輯
poly_instance =None
ring_instance =None
area_type =None
ring_temp =None
multi_type =None
group =None# 根据区域个数预先处理if';'notin poly_str:
poly_instance = ogr.Geometry(ogr.wkbPolygon)
ring_instance = ogr.Geometry(ogr.wkbLinearRing)
area_type ='single'elif poly_str =="EMPTY":continueelse:
area_type ='multiple'# 获取所有的环if( area_type =='single'):# 按照逗号分割字符串
poly_instance.AddGeometry(extract_single_ring(poly_str))elif(area_type =='multiple'):
multi_list = poly_str.split(";")
ogr_ring_list =[]for ring_list in multi_list:
ogr_ring_list.append(extract_single_ring(ring_list))# SORT THE LIST BY AREA
sorted_ogr_ring_list =sorted(ogr_ring_list, key=lambda ring: ring.GetArea(), reverse=True)
group = group_rings(sorted_ogr_ring_list)
multi_type ='hole_poly'iflen(group)==1else'multipolygon'if(multi_type =='hole_poly'):# TODO deal
poly_instance = ogr.Geometry(ogr.wkbPolygon)for item in group[0]:
poly_instance.AddGeometry(item)elif(multi_type =='multipolygon'):# TODO deal
poly_instance = ogr.Geometry(ogr.wkbMultiPolygon)
single_poly_inst = ogr.Geometry(ogr.wkbPolygon)for poly in group:for item in poly:
single_poly_inst.AddGeometry(item)
poly_instance.AddGeometry(single_poly_inst)
feature.SetGeometry(poly_instance)
layer.CreateFeature(feature)# push back the feature item
feature =None# Clear the feature to reuse it# Save the changes and close the data source
data_source =None