基于alphashape算法提取Geo点集外轮廓+ST_ChaikinSmoothing（PGSQL+PostGIS）算法对外边界平滑处理

import pickle
import psycopg2
import alphashape
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Point

# table_exist
def table_exist(table_name=None, conn=None, cur=None):
    try:
        cur.execute("select to_regclass(" + "\'" + table_name + "\'" + ") is not null")
        rows = cur.fetchall()
    except Exception as e:
        rows = []
        conn.close()
    if rows:
        data = rows
        flag = data[0][0]
        return flag
# X2GDF
def X2GDF(X, input_crs = "EPSG:2436", output_crs = "EPSG:4326", plot=True):
    # input_crs = "EPSG:2436"# input_crs = "EPSG:4326"
    # output_crs = "EPSG:4326"# output_crs = "EPSG:2436"
    coords, IDs = [], []
    for i in range(X.shape[0]):
        coords.append(Point(X[i, 0], X[i, 1]))# (X, Y)
        IDs.append("%d" % i)
    data = {'ID': IDs, 'geometry': coords}
    gdf = gpd.GeoDataFrame(data, crs=input_crs)
    # print(gdf.crs)  # 查看数据对应的投影信息(坐标参考系)
    if output_crs is not None:
        gdf.to_crs(crs=output_crs)# (Lon, Lat)

    if plot:
        gdf.plot()
        plt.show()#展示

    return gdf
# GDF2X
def GDF2X(gdata, to_epsg=2436, plot=True):
    print('input gdf crs:', gdata.crs)  # 查看数据对应的投影信息(坐标参考系) ---> epsg == 4326  (WGS_84)
    gdata = gdata.to_crs(epsg=to_epsg)# epsg == 2436 (X, Y), 4326 (lon, lat)
    print('output X crs:', gdata.crs)
    # print(gdata.columns.values)# 列名
    # print(gdata.head())  # 查看前5行数据
    if plot:
        gdata.plot()
        plt.show()#展示
    geo = gdata.geometry.values
    X0, Y0 = geo.x, geo.y
    X = np.zeros(shape=(X0.shape[0], 2), dtype=float)
    for i in range(X0.shape[0]):
        X[i, 0], X[i, 1] = X0[i],  Y0[i]
    return X
# LonLat2XY
def LonLat2XY(X, input_crs = "EPSG:4326", output_crs = "EPSG:3857", shpfile_name = "point_shp.shp"):
    # index = ['node_'+str(node) for node in range(X.shape[0])]
    # columns = ['column_'+str(col) for col in range(X.shape[1])]
    # df = pd.DataFrame(X, index=index, columns=columns)
    coords, IDs = [], []
    for i in range(X.shape[0]):
        coords.append(Point(X[i, 0], X[i, 1]))  # (X, Y)
        IDs.append("%d" % i)

    data = {'ID': IDs, 'geometry': coords}
    gdf = gpd.GeoDataFrame(data, crs=input_crs)

    print('input crs', gdf.crs)  # 查看数据对应的投影信息(坐标参考系)
    if output_crs is not None:
        gdf.to_crs(crs=output_crs)  # (Lon, Lat)

    print('output crs', gdf.crs)  # 查看数据对应的投影信息(坐标参考系)  output_crs = "EPSG:3857" or "EPSG:2436"
    if shpfile_name is not None:# shpfile = "point_shap.shp"
        gdf.to_file(filename=shpfile_name)
        print('shpfile is saved')

    geo = gdf.geometry.values
    X0, Y0 = geo.x, geo.y
    X = np.zeros(shape=(X0.shape[0], 2), dtype=float)
    for i in range(X0.shape[0]):
        X[i, 0], X[i, 1] = X0[i], Y0[i]

    return X
# get_X_fromSHP
def get_X_fromSHP(filename=r'./park/BJ_park_m.shp', epsg4326_srid=2436, plot=True):
    gdata = gpd.read_file(filename=filename, bbox=None, mask=None, rows=None)#读取磁盘上的矢量文件
    # gdata = gpd.read_file(filename=r'CUsersjeshyDocumentsArcGISJ50.gdb', layer='BOUA')#读取gdb中的矢量数据
    print('Input geodata crs:', gdata.crs)  # 查看数据对应的投影信息(坐标参考系) ---> epsg == 4326  (WGS_84)
    gdata = gdata.to_crs(epsg=epsg4326_srid)# epsg == 2436 (X, Y)
    print('Output geodata crs:', gdata.crs)
    print('Column names', gdata.columns.values)# 列名
    # print(gdata.head())  # 查看前5行数据
    if plot:
        gdata.plot()
        plt.show()#展示

    geo = gdata.geometry.values
    X0, Y0 = geo.x, geo.y
    X = np.zeros(shape=(X0.shape[0], 2), dtype=float)
    for i in range(X0.shape[0]):
        X[i, 0], X[i, 1] = X0[i],  Y0[i]

    print('n_samples:', X0.shape[0])
    output = open('X.pkl', 'wb')
    pickle.dump(X, output)
    output.close()
    print('X.pkl is saved')

    from sklearn.preprocessing import StandardScaler
    X_trans = StandardScaler().fit_transform(X)
    output = open('trans_X.pkl', 'wb')
    pickle.dump(X_trans, output)
    output.close()
    print('trans_X.pkl is saved')

    return X

if __name__ == '__main__':
    # Load the geodata
    output = open('park_X.pkl', 'rb')# park_X = get_X_fromSHP(filename=r'./park/BJ_park_m.shp', epsg4326_srid=2436, plot=True)
    X = pickle.load(output)
    output.close()

    input_crs = "EPSG:2436"
    output_crs = None

    X_uniques = np.unique(X, axis=0)# 去除重复行
    if X_uniques.shape[0] > 2:
        # get the outlines by alphashape method
        gdf_alpha_shape = alphashape.alphashape(X2GDF(X_uniques, input_crs=input_crs, output_crs=output_crs, plot=None))
        polyg = gdf_alpha_shape.geometry.values  # GeometryArray
        polyg = polyg[0]  # Not smoothed poly_line

        # pgsql-postgis with ST_ChaikinSmoothing
        Polyg_wkt_str, SRID_str = polyg.wkt, input_crs.split(':')[-1]
        sql = 'select ST_ChaikinSmoothing(ST_GeomFromEWKT(\'SRID={0};{1}\'),5) as the_geom;'.format(SRID_str, Polyg_wkt_str)

        # connection the database
        table_name = "beijing"
        conn = psycopg2.connect(database=table_name, user="postgres", password="jiangshan", host="localhost", port="5432")

        # cur = conn.cursor()
        # # CREATE TABLE IF table IS NOT EXIST
        # # 查询出来的表是否存在的状态，存在则为True，不存在则为False
        # table_flg = table_exist(table_name, conn, cur)
        # if table_flg is False:
        #     sql = "DROP TABLE public.{0} CASCADE".format(table_name)  # -- 删除表
        #     sql = "CREATE TABLE IF NOT EXISTS {0} (serial_number BIGINT, code_company TEXT, lon DOUBLE PRECISION, lat DOUBLE PRECISION, state INT)".format(table_name)
        #     cur.execute(sql)
        #     conn.commit()
        #     cur.close()

        out = gpd.read_postgis(sql=sql, con=conn, geom_col="the_geom")
        polyg = out['the_geom'].values  # GeometryArray
        polyg = polyg[0]  # Smoothed poly_line
        # close
        conn.close()