在Python中根据地理数据框内的点创建地图边界 [英] Create map boundaries from points within a geodataframe in python

问题描述

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

	几何	Cost_Trainstation_Name
1	point(1，1)	Station_1
2	点(10，10)	Station_2
...	...	...

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

这些点是使用最近邻算法创建的，因此它们不会彼此相交。

我还有一个用于整个国家的边界地理数据框架，称为boundary，我认为可能需要它来定义此地图的外部边界。和一份火车站的档案。

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

推荐答案

• 使用英国的医院作为数据来源(有效地站和医院并可互换)
• 具有英格兰边界的几何图形(以剪裁Voronoi多边形)
• 生成表示所有医院的面就变得很简单
• 使用sjoin()将面与具有所有关联属性的点相关联
• 我已使用图解来可视化和演示与之关联的面

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()

这篇关于在Python中根据地理数据框内的点创建地图边界的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！