时间戳

首页 > 解决方案 > 如何找到全球所有地区的边界框数据？

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

问题描述

我正在寻找国家内部存在的各个州的边界框数据。例如，对于印度 - 安得拉邦、卡纳塔克邦等州或挪威 - 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 中提取它吗？我正在寻找任何方式或方法来获取这些数据。即使是指出我正确方向的评论也会有所帮助。谢谢。

标签： gisbounding-box

解决方案

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

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

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

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

natural_earth_1_10_states_and_provinces

此 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://stackoverflow.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

推荐阅读