如何找到全球所有地区的边界框数据？

Posted 2023-03-27

【中文标题】如何找到全球所有地区的边界框数据？

【英文标题】：How to find bounding box data for all regions in the world?

【发布时间】：2021-11-27 09:31:32

【问题描述】：

我正在寻找存在于国家/地区内的各个州的边界框数据。例如，对于印度 - 安得拉邦、卡纳塔克邦等州或挪威 - Akershus、Aust-Agder 等的 bbox。

我知道的一种方法是从 Geofabrik 下载摘录并使用 osmium-tool 使用此命令从中获取边界框数据 -

osmium fileinfo -e -g data.bbox victoria-latest.osm.pbf

但是，Geofabrik 没有许多国家/地区的州，例如印度、中国、印度尼西亚、奥地利等。

一种替代方法是使用来自 download.openstreetmap.fr 的摘录。不幸的是，他们的提取物具有非常不正确的边界框值。例如，澳大利亚维多利亚州的边界框覆盖了澳大利亚的一半。我在他们的许多提取物中都注意到了这个问题。不丹的边界框覆盖了 3 倍的横向区域。马来西亚的提取物覆盖了印度洋的一半。这使得它们的提取无法用于从中提取边界框数据。

除了这种方式，还有这样的数据集—— https://gist.github.com/graydon/11198540

但是，如您所见，这些不包含州和地区。

我在哪里可以找到这些数据？或者，如果可能的话，我可以使用 Overpass 以某种方式从 OSM 中提取它吗？我正在寻找任何方式或方法来获取这些数据。即使是指出我正确方向的评论也会有所帮助。谢谢。

【问题讨论】：

正如我们在 gis.se (gis.stackexchange.com/questions/413305/…) 上所说的，您可以为此使用自然地球数据，也可以在 opendata.stackexchange.com 上提问

您好@IanTurton，感谢您的评论。回复您的原始评论 - 我并不是专门寻找开放数据。我正在寻找这些数据的任何来源。我已经知道的一种方法是从 Geofabrik 下载提取并使用 osmium-tool 来提取边界框数据。但是，Geofabric 提取物没有针对印度、印度尼西亚、中国等许多国家/地区的区域。不幸的是，OSM-FR 提取物的边界框数据不正确。总之，这不是一个专门针对开放数据源的问题，而是找到不同区域边界框数据的任何方法。

请下载 NaturalEarth 数据 - naturalearthdata.com/downloads/110m-cultural-vectors

谢谢。在您最初发表评论后，我开始研究它。我不熟悉 shp 文件，所以我想弄清楚如何从中提取所需的数据。

@IanTurton，我查看了 QGIS 中自然地球下载的所有文件，除了美国，这些文件中不存在其他国家/地区的州/省数据。我犯错了吗？请指导。谢谢。

【参考方案1】：

我在 Natural Earth 数据中找到了所有国家/地区/地区/省所需的数据 - https://www.naturalearthdata.com/。

但是，在使用它之后，我意识到这个数据已经很老了，对于许多国家来说，它有州但没有新的州或省（例如，印度的 Telengana 于 2014 年作为一个新州成立，但该数据集中不存在）。

1:110m 仅具有国家边界和美国各州的边界。 1:50m 有北美、澳大利亚和其他一些地区的国界 + 州界，但没有全世界。 1:10m 下载包含世界上所有国家/地区的状态。下载这个文件 - https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces.zip

要获得更新和详细的版本，请使用此摘录中的 1 级形状文件-https://biogeo.ucdavis.edu/data/gadm3.6/gadm36_levels_shp.zip

这个 python 脚本提取世界上所有国家的所有州的边界框值，并将其存储在一个 csv 文件中。

import fiona
import json
import numpy as np
# sudo pip3 install git+https://github.com/ulikoehler/UliEngineering.git
from UliEngineering.Math.Coordinates import BoundingBox
import polyline as pl
import csv

# https://gis.stackexchange.com/a/113808/138032 - how to read shape file
with fiona.open("ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp") as c:

    with open('bbox_file.csv', mode='w') as bbox_file:

        bbox_writer = csv.writer(bbox_file)
        bbox_writer.writerow(['country', 'region', 'bbox'])

        for record in c:
            country = record['properties']['admin']
            region = record['properties']['name']  # state/province name
            print(region)

            # json_object = json.dumps(record['properties'], indent=4)
            # print(json_object)
            # print(record['geometry']['type'])  # All features are either "Polygon" or "Multipolygon"

            if record['geometry']['type'] == "Polygon":
                """
                Polygon always has array(array(coordinates)). The first array is always a length of 1. So just
                accessing the array(coordinates)
                """
                coordinates = record['geometry']['coordinates'][0]
            else:
                """
                A multipolygon is basically a list of polygons, so add one more layer to it, meaning
                array(array(array(coordinates))). Here, the first layer of arrays might have any number of
                polygons inside it. And then its the same as polygons where the first array always have a 
                single element. So in array(array(array(coordinates))) - the middle array element always
                has a single element.  
                
                The current format is this -  
                [ [[(lat1, lng1), (lat2, lng2)]], [[(lat3, lng3), (lat4, lng4)]] ]
                And I need to convert it to this -  
                [(lat1, lng1), (lat2, lng2), (lat3, lng3), (lat4, lng4)]
                Just one list of coordinates. So basically, to array(coordinates). Same format as polygon. 
                
                Calling sum() for the first time, makes it into this
                [ [(lat1, lng1), (lat2, lng2)], [(lat3, lng3), (lat4, lng4)] ]
                So one layer or arrays is removed. Sum needs to be called again in order to remove another layer, 
                then it becomes this
                [(lat1, lng1), (lat2, lng2), (lat3, lng3), (lat4, lng4)]
                https://***.com/a/716489/3090120
                
                Now all the coordinates exist in a single array. 
                """
                coordinates = sum(record['geometry']['coordinates'], [])
                coordinates = sum(coordinates, [])

            # print(coordinates)  # the final single array containing all coordinates

            """
            Now the problem with these coordinates is that they are in the (longitude, latitude) format
            but I need it in the (latitude, longitude) format. So creating a new array and storing the 
            tuples in flipped format.
            """
            fixed_coordinates_array = []
            for coord_pair in coordinates:
                fixed_coordinates_array.append(coord_pair[::-1])

            # print(fixed_array)

            # https://techoverflow.net/2017/02/23/computing-bounding-box-for-a-list-of-coordinates-in-python/
            bbox = BoundingBox(np.asarray(fixed_coordinates_array))
            bbox_str = str(bbox.minx) + "," + str(bbox.miny) + "," + str(bbox.maxx) + "," + str(bbox.maxy)
            print(bbox_str)

            # coordinates to polyline - https://pypi.org/project/polyline/
            polyline = pl.encode(fixed_coordinates_array)
            # print(polyline)

            bbox_writer.writerow([country, region, bbox_str])

此外，如果您想将折线添加到 csv，只需进行这些更改 -

bbox_writer.writerow(['country', 'region', 'bbox', 'polyline'])
..
..
bbox_writer.writerow([country, region, bbox_str, polyline])

这是按国家和地区排序的最终 csv 输出 - https://gist.github.com/kuwapa/a002b7abbeeaaa1b27fd31ac9840a5dd

如果有人想更深入地为每个国家/地区的每个州/省的每个县生成边界框，他们可以使用从这里导出的形状文件 - https://gadm.org/download_world.html