在邮编级别处理地理空间数据

如何将“点”邮政编码与区域普查数据关联起来

在一些国家，邮政编码是点或路线，而不是区域。例如，加拿大邮政编码的最后三位数字对应于本地投递单位，该单位可能对应于街道一侧或农村路线上的房屋。类似地，英国邮政编码的形式为“YO8 9UR”。这可以是伦敦的单个建筑物。在5+4位的美国邮政编码中，最后四个数字确定了一个邮政投递路线（即一组地址），而不是一个区域。与常见的观点相反，美国的5位邮政编码也不是区域，它们只是一系列5+4位的邮政路线，通常从一个单独的邮局提供服务。

作为度量系统的发起者，法国非常合乎逻辑。在法国，邮政编码对应一个区域，最后两位数字对应了区（arrondissement），因此75008对应巴黎的第八区，确实是一个区域。不过，邮政投递路线可能不够理想。

由于人们和商店都有地址，并且地址与邮政编码相关联，所以大多数消费者数据都是以邮政编码为单位报告的。为了进行区域覆盖、市场份额等计算，必须确定邮政编码的区域范围。在法国这很容易，但在任何邮政编码是邮政路线而不是区域的国家将会很困难。

由于它们的邮政编码是邮政投递地址，因此可以以无限多种方式绘制多边形来划分英国/加拿大/美国为有效的邮政编码“区域”。这就是为什么英国的人口统计数据由国家统计局（ONS）按行政区域（如县）发布，而不是按邮政编码发布。美国普查数据以“邮政编码制表区”（ZCTA）级别发布，美国投票数据以县级别发布。在处理英国/加拿大/美国的数据时，您通常会有类似的地址（点）和覆盖区域收集的空间数据。如何将它们关联在一起呢？

为了说明问题，本文将结合英国邮政编码数据和普查数据进行讲解。

下载链接

如果您急于使用，您可以从https://github.com/lakshmanok/lakblogs/tree/main/uk_postcode下载我的分析结果 – 那里有几个CSV文件，它们包含您可能需要的数据。

ukpopulation.csv.gz包含以下列：

postcode,latitude,longitude,area_code,area_name,all_persons,females,males

ukpostcodes.csv.gz有一个额外的列 – 每个邮政编码的多边形以WKT格式表示：

postcode,latitude,longitude,area_code,area_name,all_persons,females,males,geometry_wkt

请注意，使用数据或代码存在风险 – 它是按“原样”分发的，没有任何明示或暗示的保证或条件。

在本文中，我将逐步介绍我如何在GitHub存储库中创建数据集。您可以使用笔记本uk_postcodes.ipynb和我一起进行操作。

原始数据

我们从三个源头的原始数据开始，这些数据根据英国开放政府许可发布：

1. Free Map Tools提供了一个可下载的文件，其中包含每个邮政编码的中心点纬度和经度。这对于空间分析来说是不够的，因为仅凭一个点无法执行诸如ST_INTERSECTS之类的操作。但它是一个很好的开始。
2. ONS已经发布了诸如“达勒姆郡”之类的区域的普查数据。但这些不是邮政编码。它们通常是在选区、郡或地区级别上。
3. 英国统计局已经帮助确定了与每个邮政编码相关联的区域。每个邮政编码属于不同的区域，这些区域为不同的目的和分辨率定义。这包括但不限于教区、选区、郡和地区。为了完整起见，这里列出了所有其他可用的关联（根据您的空间数据集，您可能需要其他列）：

pcd，pcd2，pcds，dointr，doterm，oscty，ced，oslaua，osward，parish，usertype，oseast1m，osnrth1m，osgrdind，oshlthau，nhser，ctry，rgn，streg，pcon，eer，teclec，ttwa，pct，itl，statsward，oa01，casward，park，lsoa01，msoa01，ur01ind，oac01，oa11，lsoa11，msoa11，wz11，sicbl，bua11，buasd11，ru11ind，oac11，lat，long，lep1，lep2，pfa，imd，calncv，icb，oa21，lsoa21，msoa21

我的笔记本使用wget下载数据文件：

mkdir -p indatacd indataif [ ! -f census2021.xlsx ]; then    wget -O census2021.xlsx https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/populationandhouseholdestimatesenglandandwalescensus2021/census2021/census2021firstresultsenglandwales1.xlsxfi

读取数据

直接将CSV文件读入Pandas很简单：

import pandas as pdpd.read_csv(POSTCODES_CSV)

这给我提供了每个邮政编码的中心纬度和经度：

有许多软件包可以让您将Excel文件读入Pandas，但我决定使用DuckDB，因为我稍后将在笔记本中使用它来使用SQL联接这三个数据集：

import duckdbconn = duckdb.connect()ukpop = conn.execute(f"""install spatial;load spatial;SELECT   * FROMst_read('{POPULATION_XLS}', layer='P01')""").df()

这个Excel文件有7行标题信息，我可以删除。我还将列重命名为有意义的变量：

ukpop = ukpop.drop(range(0,7), axis=0).rename(columns={    'Field1': 'area_code',     'Field2': 'area_name',     'Field3': 'all_persons',     'Field4': 'females',     'Field5': 'males'})

那是名为P01的工作表。请注意，P04工作表中有人口密度信息，但由于人口在区域代码上分布不均匀，所以这不是有用的。我们将推导出每个邮政编码的人口。

我将其写入CSV文件中，以便我可以轻松地从DuckDB中读取它。

ukpop.to_csv("temp/ukpop.csv", index=False)

类似地，我从英国统计办公室文件中提取必要的列，并将其写入CSV文件：

onspd = pd.read_csv(ONSPD_CSV)onspd = onspd[['pcd', 'oslaua', 'osward', 'parish']]onspd.to_csv("temp/onspd.csv", index=False)

关联数据

现在，我们可以使用DuckDB将这三个准备好的数据集连接起来，以获得每个邮政编码的人口密度。为什么使用DuckDB？虽然我可以在Pandas中进行连接，但我发现SQL更易读。此外，这给我一个使用新热门工具的借口。

我通过先将它们读入DuckB使用read_csv_auto来连接数据集。然后，我查找邮政编码所在的区域（parish、ward或county）的选区、教区和县，并找到报告人口密度数据的区域（parish、ward或county）：

/* pcd,oslaua,osward,parish */WITH onspd AS (    SELECT       *     FROM    read_csv_auto('temp/onspd.csv', header=True)),/* area_code,area_name,all_persons,females,males */ukpop AS (    SELECT       *     FROM    read_csv_auto('temp/ukpop.csv', header=True)),/* id,postcode,latitude,longitude */postcodes AS (    SELECT       *     FROM    read_csv_auto('indata/ukpostcodes.csv', header=True)),/* postcode, area_code */postcode_to_areacode AS (  SELECT     pcd AS postcode,    ANY_VALUE(area_code) as area_code  FROM onspd  JOIN ukpop   ON (area_code = oslaua OR area_code = osward OR area_code = parish)  GROUP BY pcd)SELECT  postcode, latitude, longitude, /* from postcodes */  area_code, area_name, /* from ukpop */  all_persons,females,males /* from ukpop, but has to be spatially corrected */FROM postcode_to_areacodeJOIN postcodes USING (postcode)JOIN ukpop USING (area_code)

请注意，空间数量是指整个区域而不是邮政编码的标量。它们必须在邮政编码之间进行划分。

在邮政编码之间划分区域数量

all_persons、females、males都对应整个区域，而不是特定的邮政编码。我们可以根据邮政编码的区域面积按比例划分，但可以适应邮政编码的无限多个多边形，而且后面我们将看到，靠近公园和湖泊的邮政编码的面积范围有些模糊。因此，我们将采取一些简单的方法，以获得一个唯一的答案-我们将在整个区域的所有邮政编码之间均匀分配标量值！这听起来并不奇怪-在密度较高的社区中，有更多的邮政编码，因此在邮政编码之间进行均等划分相当于按照人口密度比例分配标量数量。

npostcodes = ukpop.groupby('area_code')['postcode'].count()for col in ['females', 'males', 'all_persons']:    ukpop[col] = ukpop.apply(lambda row:  row[col]/npostcodes[row['area_code']], axis=1)

到这一步，我们已经得到了每个邮政编码的数量-这就是我们需要的关联：

所以，写出来：

ukpop.to_csv("ukpopulation.csv", index=False)

邮政编码的面积范围

对于许多分析，我们希望邮政编码不是点而是区域。尽管我们可以使用无限多个多边形来划分英国，以使每个多边形中只有一个邮政编码中心点，但确实存在一个“最佳”多边形的概念。那就是Voronoi分区，它将区域划分为使任何点都属于最近的邮政编码：

为了计算这个，我们可以使用scipy：

import numpy as npfrom scipy.spatial import Voronoi, voronoi_plot_2dpoints = df[['latitude', 'longitude']].to_numpy()vor = Voronoi(points)

我在这里假设区域足够小，以至于从纬度和经度计算出的大地距离和欧几里得距离之间没有太大的区别。英国的邮政编码足够小，这种情况是成立的。

结果的组织方式是对于每个点，都有一个由一组顶点组成的区域。我们可以使用以下代码为每个点创建一个WKT多边形字符串：

def make_polygon(point_no):    region_no = vor.point_region[point_no]    region = vor.regions[region_no]    if len(region) >= 3:        # close the ring        closed_region = region.copy()        closed_region.append(closed_region[0])        # create a WKT of the points        polygon = "POLYGON ((" + ','.join([ f"{vor.vertices[v][1]} {vor.vertices[v][0]}" for v in closed_region]) + "))"        return polygon    else:        return None

下面是一个示例结果：

POLYGON ((-0.32491691953979235 51.7393550489536,-0.32527234008402217 51.73948967705648,-0.32515738641624575 51.73987124225542,-0.3241646650618929 51.74087626616231,-0.3215663358407994 51.742660660928614,-0.32145633473723817 51.742228570262824,-0.32491691953979235 51.7393550489536))

我们可以创建一个GeoDataFrame并绘制一部分邮政编码：

import geopandas as gpdfrom shapely import wktdf['geometry'] = gpd.GeoSeries.from_wkt(df['geometry_wkt'])gdf = gpd.GeoDataFrame(df, geometry='geometry')gdf[gdf['area_name'] == 'St Albans'].plot()

这是伯明翰：

gdf[gdf['area_name'] == 'Birmingham'].plot()

无人居住区域

请注意顶部的角和中间的大片蓝色区域。发生了什么？让我们看看 Google 地图上的伯明翰：

注意公园区域吗？皇家邮政不必向那里的任何人投递邮件。因此，那里没有邮政编码。因此，附近的邮政编码会“扩展”到这些区域中。这将导致空间计算中的问题，因为这些邮政编码看起来比它们实际上大得多。

为了解决这个问题，我将采用一种相当启发式的方法。我将将英国划分为0.01×0.01（约1平方公里）分辨率的网格单元，并找到其中没有邮政编码的网格单元：

GRIDRES = 0.01min_lat, max_lat = np.round(min(df['latitude']), 2) - GRIDRES, max(df['latitude']) + GRIDRESmin_lon, max_lon = np.round(min(df['longitude']), 2) - GRIDRES, max(df['longitude']) + GRIDRESprint(min_lat, max_lat, min_lon, max_lon)npostcodes = np.zeros([ int(1+(max_lat-min_lat)/GRIDRES), int(1+(max_lon-min_lon)/GRIDRES) ])for point in points:    latno = int((point[0] - min_lat)/GRIDRES)    lonno = int((point[1] - min_lon)/GRIDRES)    npostcodes[latno, lonno] += 1unpop = []for latno in range(len(npostcodes)):    for lonno in range(len(npostcodes[latno])):        if npostcodes[latno][lonno] == 0:            # 没有人住在这里。            # 为这个位置编写一个邮政编码            # 邮编 纬度 经度 区域代码 区域名称 每平方公里人数            unpop.append({                'postcode': f'UNPOP {latno}x{lonno}',                'latitude': min_lat + latno * 0.01,                'longitude': min_lon + lonno * 0.01,                'all_persons': 0            })               

我们将在这些无人居住的网格单元中创建假的邮政编码，并将人口密度设置为零。将这些假的邮政编码添加到实际的邮政编码中，并重复 Voronoi 分析：

df2 = pd.concat([df, pd.DataFrame.from_records(unpop)])points = df2[['latitude', 'longitude']].to_numpy()vor = Voronoi(points)df2['geometry_wkt'] = [make_polygon(x) for x in range(len(vor.point_region))]df2['geometry'] = gpd.GeoSeries.from_wkt(df2['geometry_wkt'])gdf = gpd.GeoDataFrame(df2, geometry='geometry')

现在，当我们绘制伯明翰时，我们得到了一个更好的结果：

这就是我将作为第二个 CSV 文件保存的数据帧：

gdf.to_csv("ukpostcodes.csv", index=False)

[可选] 加载到 BigQuery

我们可以将 CSV 文件加载到 BigQuery 中，并对其进行一些空间分析，但最好先让 BigQuery 将最后一列字符串解析为几何图形，并按邮政编码对数据进行聚类：

CREATE OR REPLACE TABLE uk_public_data.postcode_popgeo2
CLUSTER BY postcode
AS
SELECT * EXCEPT(geometry_wkt),
  SAFE.ST_GEOGFROMTEXT(geometry_wkt, make_valid=>TRUE) AS geometry,
FROM uk_public_data.postcode_popgeo

现在，我们可以轻松查询它。例如，我们可以使用 ST_AREA 获取邮政编码的面积：

SELECT COUNT(*) AS num_postcodes,
  SUM(ST_AREA(geometry))/1e6 AS total_area,
  SUM(all_persons) AS population,
  area_name
FROM uk_public_data.postcode_popgeo2
GROUP BY area_name
ORDER BY population DESC

摘要

空间分析通常需要面积范围，而不仅仅是点位置。在邮政编码是点/线路的国家，可以通过无限多种方式生成邮政编码的多边形空间范围。一个合理的方法是使用 Voronoi 区域创建包含这些邮政编码的多边形。然而，如果这样做，你将在湖泊或公园附近得到不自然地大的多边形，邮局不会投递邮件。为了解决这个问题，还需要对国家进行网格划分，并在未有人居住的网格单元中创建人工邮政编码。在本文中，我演示了如何为英国进行这样的操作。相关的笔记本可以用于其他地方。

下一步

1. 在 https://github.com/lakshmanok/lakblogs/blob/main/uk_postcode/uk_postcodes.ipynb 上查看完整的代码
2. 从 https://github.com/lakshmanok/lakblogs/tree/main/uk_postcode 下载数据