from osgeo import ogr
import pandas as pd
def extract_single_ring( point_list_str ):
ring_instance = ogr.Geometry(ogr.wkbLinearRing)
point_list = point_list_str.split(',')
count = 0
first_x = None
first_y = None
for 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.
"""
def group_rings(ogr_ring_list):
grouped_rings = []
while(ogr_ring_list):
current_ring = ogr_ring_list.pop(0)
group = [current_ring]
i = 0
while 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 ';' not in poly_str:
poly_instance = ogr.Geometry(ogr.wkbPolygon)
ring_instance = ogr.Geometry(ogr.wkbLinearRing)
area_type = 'single'
elif poly_str == "EMPTY":
continue
else:
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' if len(group) == 1 else '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
csv2geojson
最新推荐文章于 2025-01-09 09:58:48 发布
本文介绍了一个Python函数,用于根据面积对地理空间中的线性环进行分组,最后将结果转换为GeoJSON格式。函数首先处理单个环,然后处理多个环并按面积排序,形成多边形或包含洞的多边形。
400

被折叠的 条评论
为什么被折叠?



