如何在特定区域为 H3 六边形生成 shapefile
Posted
技术标签:
【中文标题】如何在特定区域为 H3 六边形生成 shapefile【英文标题】:How to generate shapefiles for H3 hexagons in a particular area 【发布时间】:2018-12-12 01:02:17 【问题描述】:我想为特定地理区域中的H3 六边形生成 shapefile。特别是,我对分辨率为 6、7 和 9 的湾区感兴趣。如何为覆盖该区域的六边形创建 shapefile?
我是 shapefile 或任何其他地理数据结构的新手。我最喜欢 python 和 R。
【问题讨论】:
【参考方案1】:这里的基本步骤是:
选取所需区域的多边形。边界框应该可以正常工作。 使用polyfill
方法以所需分辨率用六边形填充多边形。
遍历每个六边形并使用h3ToGeoBoundary
函数获取边界。
将这些边界放入 GeoJSON 文件中
使用ogr2ogr
之类的转换器转换为shapefile。
Python 绑定尚未发布,我对 R 绑定不熟悉,但 javascript 版本可能如下所示:
var h3 = require('h3-js');
var bbox = [
[-123.308821530582, 38.28055644998254],
[-121.30037257250085, 38.28055644998254],
[-121.30037257250085, 37.242722073589164],
[-123.308821530582, 37.242722073589164]
];
var hexagons = h3.polyfill(bbox, 6, true);
var geojson =
type: 'Feature',
geometry:
type: 'MultiPolygon',
coordinates: hexagons.map(function toBoundary(hex)
return [h3.h3ToGeoBoundary(hex, true)];
)
;
console.log(JSON.stringify(geojson));
你会像这样使用脚本:
node bbox-geojson.js | ogr2ogr -f "ESRI Shapefile" bbox-hexagons.shp /vsistdin/
【讨论】:
不错的脚本。谢谢。您认为有什么方法可以将索引保留为每个多边形的属性? @loki - 我相信你需要创建一个FeatureCollection
的 Polygon
功能,每个功能都可以有一个 id
和可能有一个 properties
【参考方案2】:
如果您在R
中寻找解决方案,h3jsr
package 提供对 Uber 的 H3 库的访问。可以使用函数h3jsr::polyfill()
和h3jsr::h3_to_polygon
来解决您的问题。
可重现的例子
library(ggplot2)
library(h3jsr)
library(sf)
library(sf)
# read the shapefile of the polygon area you're interested in
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
# projection
nc <- st_transform(nc, crs = 4326)
# get the unique h3 ids of the hexagons intersecting your polygon at a given resolution
nc_5 <- polyfill(nc, res = 5, simple = FALSE)
# pass the h3 ids to return the hexagonal grid
hex_grid5 <- unlist(nc_5$h3_polyfillers) %>% h3_to_polygon(simple = FALSE)
这将返回下面的多边形:
【讨论】:
你有机会把它变成 Python 的 H3 吗?截至今天,我无法使用 R 进行复制。【参考方案3】:在这里接受 John Stud 的问题,因为我遇到了同样的“问题”。在下文中,我将评论如何读取 shapefile,使用 H3 对其进行六边形化,并从中获取 Hexagon 地理数据框(并最终将其保存为 shapefile)。
可重现的示例
让我们为美国获取一个 shapefile,例如here(我用的是“cb_2018_us_state_500k.zip”那个)。
# Imports
import h3
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import shapely
from shapely.ops import unary_union
from shapely.geometry import mapping, Polygon
# Read shapefile
gdf = gpd.read_file("data/cb_2018_us_state_500k.shp")
# Get US without territories / Alaska + Hawaii
us = gdf[~gdf.NAME.isin(["Hawaii", "Alaska", "American Samoa",
"United States Virgin Islands", "Guam",
"Commonwealth of the Northern Mariana Islands",
"Puerto Rico"])]
# Plot it
fig, ax = plt.subplots(1,1)
us.plot(ax=ax)
plt.show()
# Convert to EPSG 4326 for compatibility with H3 Hexagons
us = us.to_crs(epsg=4326)
# Get union of the shape (whole US)
union_poly = unary_union(us.geometry)
# Find the hexagons within the shape boundary using PolyFill
hex_list=[]
for n,g in enumerate(union_poly):
if (n+1) % 100 == 0:
print(str(n+1)+"/"+str(len(union_poly)))
temp = mapping(g)
temp['coordinates']=[[[j[1],j[0]] for j in i] for i in temp['coordinates']]
hex_list.extend(h3.polyfill(temp,res=5))
# Create hexagon data frame
us_hex = pd.DataFrame(hex_list,columns=["hex_id"])
# Create hexagon geometry and GeoDataFrame
us_hex['geometry'] = [Polygon(h3.h3_to_geo_boundary(x, geo_json=True)) for x in us_hex["hex_id"]]
us_hex = gpd.GeoDataFrame(us_hex)
# Plot the thing
fig, ax = plt.subplots(1,1)
us_hex.plot(ax=ax, cmap="prism")
plt.show()
上图分辨率为“5”(https://h3geo.org/docs/core-library/restable/),建议你也看看其他分辨率,比如4:
当然,这取决于“缩放级别”,即您是查看整个国家还是仅查看城市等。
当然,要回答最初的问题:您可以使用
保存生成的 shapefileus_hex.to_file("us_hex.shp")
【讨论】:
【参考方案4】:修改了@nrabinowitz 中的一个,带有单独的多边形和索引名称:
var h3 = require('h3-js');
var bbox = [
[-123.308821530582, 38.28055644998254],
[-121.30037257250085, 38.28055644998254],
[-121.30037257250085, 37.242722073589164],
[-123.308821530582, 37.242722073589164]
];
var hexagons = h3.polyfill(bbox, 5, false);
var features = hexagons.map(function toBoundary(hex)
var coords = h3.h3ToGeoBoundary(hex, true)
var feature = "type": "Feature",
"properties": "name": hex,
"geometry":
"type": "Polygon",
"coordinates": [coords];
return feature;
);
console.log(JSON.stringify(
"type": "FeatureCollection",
"features": features));
【讨论】:
以上是关于如何在特定区域为 H3 六边形生成 shapefile的主要内容,如果未能解决你的问题,请参考以下文章
随机生成 2d 多边形/区域以填充 [长度,宽度] 区域的方法?