csv2geojson

本文介绍了一个Python函数,用于根据面积对地理空间中的线性环进行分组,最后将结果转换为GeoJSON格式。函数首先处理单个环,然后处理多个环并按面积排序,形成多边形或包含洞的多边形。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值