从python中的地理数据框中的点创建地图边界

Question

我有一个名为map 的地理数据框，其中包含一个点列表和一个列Closest_TrainStation_name，其中包含离该点最近的火车站的名称。

	geometry	Closest_TrainStation_name
1	POINT(1,1)	Station_1
2	POINT(10,10)	Station_2
...	...	...

是否可以创建包含每个组的边界多边形，如下图所示？每个多边形都有离原始文件最近的火车站的名称。

这些点是使用最近邻算法制作的，因此它们不会相互交叉。

我还有一个名为boundary 的全国边界地理数据框，我认为可能需要它来定义此地图的外部边界。还有火车站的档案。

我能找到的所有方法都是关于从边界上的点创建边界，例如凸包。

Answer 1

• 已将英格兰的医院用作数据源（有效地车站和医院并且可互换）
• 具有英格兰边界的几何图形（以裁剪 Voronoi 多边形）
• 然后生成代表所有医院的多边形​​变得简单
• 使用sjoin() 将多边形与带有所有相关属性的点相关联
• 已使用 plotly 来可视化和演示具有与其相关联的医院属性的多边形

import shapely.ops

# points for hospitals
gdf = gpd.GeoDataFrame(
    geometry=dfhos.loc[:, ["Longitude", "Latitude"]].apply(
        shapely.geometry.Point, axis=1
    )
)

# generate voroni polygon for each hospital
gdfv = gpd.GeoDataFrame(
    geometry=[
        p.intersection(uk)
        for p in shapely.ops.voronoi_diagram(
            shapely.geometry.MultiPoint(gdf["geometry"].values)
        ).geoms
    ]
)

# spatial join polygons to points to pick up full details of hospital
gdf3 = gpd.sjoin(gdfv, gdf, how="left").merge(
    dfhos, left_on="index_right", right_index=True
)
gdf3["Color"] = pd.factorize(gdf3["Postcode"], sort=True)[0]

# and visualize
fig = (
    px.choropleth_mapbox(
        gdf3,
        geojson=gdf3.__geo_interface__,
        locations=gdf3.index,
        hover_data=["OrganisationCode","OrganisationName","Postcode"],
        color="Color",
        color_continuous_scale="phase",
    )
    .update_layout(
        mapbox={
            "style": "carto-positron",
            "center": {
                "lon": sum(gdf3.total_bounds[[0, 2]]) / 2,
                "lat": sum(gdf3.total_bounds[[1, 3]]) / 2,
            },
            "zoom": 5,
        },
        margin={"l": 0, "r": 0, "t": 0, "b": 0},
        coloraxis={"showscale":False}
    )
)

fig

获取英格兰医院的位置和英格兰边界的多边形

import geopandas as gpd
import shapely.geometry
import numpy as np
import plotly.express as px
import requests, io
from pathlib import Path
from zipfile import ZipFile
import urllib
import pandas as pd

# fmt: off
# uk geometry
url = "http://geoportal1-ons.opendata.arcgis.com/datasets/687f346f5023410ba86615655ff33ca9_1.zip"
f = Path.cwd().joinpath(urllib.parse.urlparse(url).path.split("/")[-1])

if not f.exists():
    r = requests.get(url, stream=True, headers={"User-Agent": "XY"})
    with open(f, "wb") as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)
    zfile = ZipFile(f)
    zfile.extractall(f.stem)

f2 = Path.cwd().joinpath("uk.geojson")
if not f2.exists():
    gdf2 = gpd.read_file(list(f.parent.joinpath(f.stem).glob("*.shp"))[0])
    gdf2 = gdf2.loc[gdf2["ctyua16cd"].str[0] == "E"]
    uk = gpd.GeoDataFrame(geometry=[p for p in shapely.ops.unary_union(gdf2.to_crs(gdf2.estimate_utm_crs())["geometry"].values).simplify(5000).geoms]).set_crs(gdf2.estimate_utm_crs()).to_crs("EPSG:4326")

    uk.to_file(Path.cwd().joinpath("uk.geojson"), driver='GeoJSON')
uk = gpd.read_file(f2)
uk = shapely.geometry.MultiPolygon(uk["geometry"].values)
# fmt: on

# get hospitals in UK
dfhos = pd.read_csv(io.StringIO(requests.get("https://assets.nhs.uk/data/foi/Hospital.csv").text),sep="Č",engine="python",)
dfhos = dfhos.loc[lambda d: d["Sector"].eq("NHS Sector") & d["SubType"].eq("Hospital")].groupby("ParentODSCode").first()

Answer 2

• 已经生成了一些样本数据
• 使用https://geopandas.org/en/stable/docs/user_guide/aggregation_with_dissolve.html 的组合生成MultiPoint shapely 对象，然后convex_hull 将为您提供关联点的边界多边形

import io
import pandas as pd
import geopandas as gpd
import shapely.wkt

df = pd.read_csv(
    io.StringIO(
        """geometry,stationname
POINT (6.84365570015329 53.29316918283323),station 0
POINT (6.84397077786025 53.34219781242178),station 0
POINT (6.86411903336562 53.31761406936343),station 0
POINT (6.818510924831408 53.31265038241855),station 0
POINT (6.56484740742038 53.29280081541883),station 1
POINT (6.626702808494024 53.28117175507543),station 1
POINT (6.593563692631294 53.3071306760684),station 1
POINT (6.578341784305515 53.33493210605594),station 1
POINT (6.622799032233983 53.28010824623466),station 1
POINT (6.613683981274606 53.29073387043917),station 1
POINT (7.132754969933528 53.04549586270261),station 2
POINT (7.198050554432563 53.1025440484535),station 2
POINT (7.129751539288898 53.05572095828281),station 2
POINT (7.176645494114577 53.11826239684777),station 2
POINT (7.197200682695362 53.07260800358519),station 2
POINT (7.162805175737367 53.10861252890756),station 2
POINT (7.109351059746762 53.13349222863155),station 2
POINT (7.06391508890864 53.11204313452662),station 2
POINT (7.184832998096279 53.05647939687099),station 2
POINT (7.159339780758176 53.06749816133893),station 2
POINT (7.107083795756558 53.08269279330896),station 2
POINT (7.023547595125693 53.07391566675645),station 2
POINT (7.150033380903857 53.1339122880752),station 2
POINT (7.087567734179277 53.13629412269268),station 2
POINT (7.082082576050039 53.09928180316115),station 2
POINT (7.163037435775678 53.10643603535302),station 2
POINT (7.188994818425656 53.11808478608855),station 2
POINT (7.096403367400944 53.09224186363327),station 2
POINT (7.166715752727275 53.13534692250609),station 2
POINT (7.153873895859038 53.15187182174655),station 2"""
    )
)

gdf = gpd.GeoDataFrame(df, geometry=df["geometry"].apply(shapely.wkt.loads)).dissolve(
    "stationname"
)["geometry"].apply(lambda mp: mp.convex_hull)

gdf.plot()

stationname	geometry
station 0	POLYGON ((6.84365570015329 53.29316918283323, 6.818510924831408 53.31265038241855, 6.84397077786025 53.34219781242178, 6.86411903336562 53.31761406936343, 6.84365570015329 53.29316918283323))
station 1	POLYGON ((6.622799032233983 53.28010824623466, 6.56484740742038 53.29280081541883, 6.578341784305515 53.33493210605594, 6.626702808494024 53.28117175507543, 6.622799032233983 53.28010824623466))
station 2	POLYGON ((7.132754969933528 53.04549586270261, 7.023547595125693 53.07391566675645, 7.087567734179277 53.13629412269268, 7.153873895859038 53.15187182174655, 7.188994818425656 53.11808478608855, 7.198050554432563 53.1025440484535, 7.197200682695362 53.07260800358519, 7.184832998096279 53.05647939687099, 7.132754969933528 53.04549586270261))