使用地理空间数据

本教程的这一部分讨论如何使用 geopandasshapely 在Python中操作地理空间数据。如果您以前从未使用过这些库,或者正在寻找有关它们如何工作的复习资料,则此页适合您!

我建议在本教程中使用 Binder .

坐标参考系

这个 GeoDataFramepandas DataFrame 使用附着的几何图形:

import pandas as pd; pd.set_option('max_columns', 6)  # Unclutter display.
import geopandas as gpd
import geoplot as gplt

# load the example data
nyc_boroughs = gpd.read_file(gplt.datasets.get_path('nyc_boroughs'))
nyc_boroughs
BoroCode BoroName Shape_Leng Shape_Area geometry
0 5 Staten Island 330385.03697 1.623853e+09 (POLYGON ((-74.05050806403247 40.5664220341608...
1 4 Queens 861038.47930 3.049947e+09 (POLYGON ((-73.83668274106708 40.5949466970158...
2 3 Brooklyn 726568.94634 1.959432e+09 (POLYGON ((-73.8670614947212 40.58208797679338...
3 1 Manhattan 358532.95642 6.364422e+08 (POLYGON ((-74.01092841268033 40.6844914725429...
4 2 Bronx 464517.89055 1.186804e+09 (POLYGON ((-73.89680883223775 40.7958084451597...

大多数操作将在 pandas DataFrame 将在一个 GeoDataFrame ,而处理当前几何体的方法并不多。其中最明显的是添加了一个用于存储几何图形的列,可以使用 geometry 属性:

nyc_boroughs.geometry
0    (POLYGON ((-74.05050806403247 40.5664220341608...
1    (POLYGON ((-73.83668274106708 40.5949466970158...
2    (POLYGON ((-73.8670614947212 40.58208797679338...
3    (POLYGON ((-74.01092841268033 40.6844914725429...
4    (POLYGON ((-73.89680883223775 40.7958084451597...
Name: geometry, dtype: object

无论何时,在 GeoDataFrame ,首先要检查它的 坐标参考系 .

A coordinate reference system 或CRS,是一种定义空间中点的位置的系统。您可以使用 crs 属性:

nyc_boroughs.crs
{'init': 'epsg:4326'}

In this case epsg:4326 is the official identifier for what the rest of us more commonly refer to as “longitude and latitude”. Most coordinate reference systems have a well-defined EPSG number, which you can look up using the handy spatialreference.org website.

为什么还有经纬度以外的坐标系?例如,美国地理位置服务局维护美国极高精度的地图,为美国各地区维护110个坐标参考系,称为“州平面坐标系”。经纬度使用 spherical coordinates ;州平面坐标系使用“平地” Cartesian coordinate . 因此,状态平面坐标在计算上要简单得多,同时对大多数应用程序保持足够精确(在其“区域”内)。

由于这个原因,州平面坐标系统仍然在整个政府中使用。例如,以下是从纽约市发布的MapPLUTO数据集中获取的数据样本:

nyc_map_pluto_sample = gpd.read_file(gplt.datasets.get_path('nyc_map_pluto_sample'))
nyc_map_pluto_sample
Borough Block Lot ... Shape_Leng Shape_Area geometry
0 MN 1 10 ... 12277.824113 7.550340e+06 POLYGON ((979561.8712409735 191884.2491553128,...
1 MN 1 101 ... 3940.840373 5.018974e+05 POLYGON ((972382.8255597204 190647.2667211443,...
2 MN 1 101 ... 3940.840373 5.018974e+05 POLYGON ((972428.8290766329 190679.1751218885,...
3 MN 1 101 ... 3940.840373 5.018974e+05 POLYGON ((972058.3399882168 190689.2800885588,...
4 MN 1 201 ... 6306.268341 1.148539e+06 POLYGON ((973154.7118112147 194614.3312935531,...
5 MN 2 1 ... 2721.060649 1.008250e+05 POLYGON ((980915.0020648837 194319.1402828991,...
6 MN 2 2 ... 2411.869687 8.724423e+04 POLYGON ((981169.004181549 194678.8213220537, ...

7 rows × 90 columns

这些数据存储在长岛州平面坐标系中 (EPSG 2263 _). 不幸的是,CRS on read被错误地设置为 epsg:4326 我们必须自己把它设置到正确的坐标系。

nyc_map_pluto_sample.crs = {'init': 'epsg:2263'}
nyc_map_pluto_sample.crs
{'init': 'epsg:2263'}

根据数据集, crs may be set to either epsg:<INT> or to a raw proj4 投影字典。底线是,读入数据集后,始终验证数据集坐标参照系是否设置为其文档应设置的坐标系。

如果你确定你的坐标不是经纬度,通常你要做的第一件事就是转换成它。 to_crs 这是否:

nyc_map_pluto_sample = nyc_map_pluto_sample.to_crs(epsg=4326)
nyc_map_pluto_sample
Borough Block Lot ... Shape_Leng Shape_Area geometry
0 MN 1 10 ... 12277.824113 7.550340e+06 POLYGON ((-74.0169058260488 40.69335342975063,...
1 MN 1 101 ... 3940.840373 5.018974e+05 POLYGON ((-74.04279194703045 40.68995148413111...
2 MN 1 101 ... 3940.840373 5.018974e+05 POLYGON ((-74.04262611856618 40.69003912689961...
3 MN 1 101 ... 3940.840373 5.018974e+05 POLYGON ((-74.04396208819837 40.69006636010664...
4 MN 1 201 ... 6306.268341 1.148539e+06 POLYGON ((-74.04001513069795 40.7008411559464,...
5 MN 2 1 ... 2721.060649 1.008250e+05 POLYGON ((-74.01202751677701 40.70003725302833...
6 MN 2 2 ... 2411.869687 8.724423e+04 POLYGON ((-74.01111163437271 40.70102458543801...

7 rows × 90 columns

协调顺序

shapely ,类库 geopandas 用于存储其几何图形,使用“现代”经纬度 (x, y) 协调秩序。这与“历史”经纬度不同 (y, x) 协调秩序。“在野外”的数据集可以是任意一种格式。

我们没有办法 geopandas 了解数据集在加载时是采用一种格式还是另一种格式。将数据集转换为正确的坐标系后,请务必再次检查几何图形是否也处于正确的坐标顺序。

这是一个很容易犯的错误,人们总是在犯这个错误!

确保坐标顺序正确的最快方法是知道数据的正确x坐标和y坐标应该是什么,然后仔细观察。

几何形状类型

世界的每一个元素 geometry column in a GeoDataFrame is a shapely object. Shapely 是一个几何操作库,用于在空间中操纵几何图形,是处理形状数据的首选pythonapi。

shapely 仅定义少数几种几何图形类型:

  • Point -一个要点。

  • MultiPoint -一组点。

  • LineString -线段。

  • MultiLineString -一组线(如一系列相连的线段)。

  • LinearRing -封闭的线条集合。基本上是一个面积为零的多边形。

  • Polygon -沿着一系列点的闭合形状。

  • MultiPolygon -多边形的集合。

你可以查一下 type 使用 type 操作员:

type(nyc_boroughs.geometry.iloc[0])
shapely.geometry.multipolygon.MultiPolygon
type(nyc_map_pluto_sample.geometry.iloc[0])
shapely.geometry.polygon.Polygon

执行几何运算

这个 shapely user manual 提供可使用库执行的大量几何操作列表:从简单操作(如平移和变换)到更复杂的操作(如多边形缓冲)。

您可以使用本机方法对对象中的几何体逐对象应用变换 pandas map function on the geometry column. For example, here is one way of deconstructing a set of Polygon or MultiPolygon objects into simplified convex hulls

%time gplt.polyplot(nyc_boroughs.geometry.map(lambda shp: shp.convex_hull))
CPU times: user 62.7 ms, sys: 2.64 ms, total: 65.3 ms
Wall time: 71.6 ms
<matplotlib.axes._subplots.AxesSubplot at 0x11c61d7b8>
../_images/Working_with_Geospatial_Data_18_2.png

您可以用这种方法对图形执行任意复杂的几何变换。然而, most common operations 以优化的形式提供,作为 geopandas 应用程序编程接口。下面是一种创建凸面外壳的快速方法,例如:

%time nyc_boroughs.convex_hull
CPU times: user 55.6 ms, sys: 6.45 ms, total: 62.1 ms
Wall time: 39.9 ms
0    POLYGON ((-74.24712436215984 40.49611539517034...
1    POLYGON ((-73.94073681665428 40.54182008715522...
2    POLYGON ((-73.98336058039274 40.56952999448672...
3    POLYGON ((-74.02305574749596 40.68291694544512...
4    POLYGON ((-73.87830680057651 40.78535662050845...
dtype: object

深入研究地理空间数据转换超出了本简短指南的范围。我只想说它们有很多,你可以通过咨询 geopandasshapely 文档。

定义自己的几何图形

在本教程的这一节中,我们将重点介绍 shapely 很可能会出现:定义自己的几何图形。

在上面的例子中,我们直接从地理空间文件中读取了一个GeoDataFrame:我们的自治区信息存储在 GeoJSON 格式,而我们的建筑足迹是 Shapefile . 如果我们把地理空间数据嵌入到一个普通的 CSVJSON 文件,读入一个普通的 pandas DataFrame 是吗?

nyc_collisions_sample = pd.read_csv(gplt.datasets.get_path('nyc_collisions_sample'))
nyc_collisions_sample
LATITUDE LONGITUDE DATE TIME
0 40.767373 -73.950057 04/16/2016 4:13
1 40.862670 -73.909039 04/16/2016 4:30
2 40.716507 -73.961275 04/16/2016 4:30
3 40.749788 -73.987768 04/16/2016 4:30
4 40.702401 73.960496 04/16/2016 4:50

以非地理空间格式保存包含轻型地理空间数据的数据集(例如点,可能是线段,但通常不是整个多边形)是非常常见的。

在这种情况下,您可以导入 shapely 直接使用它来定义我们自己的几何图形,然后初始化 GeoDataFrame . 这个 pandas apply 功能是最好的:

from shapely.geometry import Point

collision_points = nyc_collisions_sample.apply(
    lambda srs: Point(float(srs['LONGITUDE']), float(srs['LATITUDE'])),
    axis='columns'
)
collision_points
0           POINT (-73.950057 40.767373)
1    POINT (-73.90903900000001 40.86267)
2           POINT (-73.961275 40.716507)
3           POINT (-73.987768 40.749788)
4    POINT (73.96049599999999 40.702401)
dtype: object

从那里我们把这个几何级数传递给 geometry 性质 GeoDataFrame initializer:

import geopandas as gpd
nyc_collisions_sample_geocoded = gpd.GeoDataFrame(nyc_collisions_sample, geometry=collision_points)
nyc_collisions_sample_geocoded
LATITUDE LONGITUDE DATE TIME geometry
0 40.767373 -73.950057 04/16/2016 4:13 POINT (-73.950057 40.767373)
1 40.862670 -73.909039 04/16/2016 4:30 POINT (-73.90903900000001 40.86267)
2 40.716507 -73.961275 04/16/2016 4:30 POINT (-73.961275 40.716507)
3 40.749788 -73.987768 04/16/2016 4:30 POINT (-73.987768 40.749788)
4 40.702401 73.960496 04/16/2016 4:50 POINT (73.96049599999999 40.702401)

在大多数情况下,CSV中提供的具有地理空间信息的数据将是与单个坐标相对应的点数据。但是,有时可能希望定义更复杂的几何图形:例如,正方形区域和 也许 吧 甚至是复杂的多边形。虽然我们不讨论这些情况,但它们非常类似于我们在这里展示的非常简单的点情况。有关此任务的进一步参考,请参阅 shapely 文档。

连接现有几何图形

有时必要的地理空间数据完全在别处。

假设现在我们有了各州关于肥胖的信息。

obesity = pd.read_csv(gplt.datasets.get_path('obesity_by_state'), sep='\t')
obesity.head()
State Percent
0 Alabama 32.4
1 Missouri 30.4
2 Alaska 28.4
3 Montana 24.6
4 Arizona 26.8

我们想把这些信息放在地图上。但我们没有任何几何学!

我们将再次定义一个几何体。但这次,我们不需要编写自己的数据,而是需要找到具有状态形状的数据,并将该数据与该数据关联起来。在其他情况下,可能有其他形状:警察辖区、调查区等。以下就是这样一个数据集:

contiguous_usa = gpd.read_file(gplt.datasets.get_path('contiguous_usa'))
contiguous_usa.head()
state adm1_code population geometry
0 Minnesota USA-3514 5303925 POLYGON ((-89.59940899999999 48.010274, -89.48...
1 Montana USA-3515 989415 POLYGON ((-111.194189 44.561156, -111.291548 4...
2 North Dakota USA-3516 672591 POLYGON ((-96.601359 46.351357, -96.5389080000...
3 Idaho USA-3518 1567582 POLYGON ((-111.049728 44.488163, -111.050245 4...
4 Washington USA-3519 6724540 POLYGON ((-116.998073 46.33017, -116.906528 46...

简单的 join 解决问题:

result = contiguous_usa.set_index('state').join(obesity.set_index('State'))
result.head()
adm1_code population geometry Percent
state
Minnesota USA-3514 5303925 POLYGON ((-89.59940899999999 48.010274, -89.48... 25.5
Montana USA-3515 989415 POLYGON ((-111.194189 44.561156, -111.291548 4... 24.6
North Dakota USA-3516 672591 POLYGON ((-96.601359 46.351357, -96.5389080000... 31.0
Idaho USA-3518 1567582 POLYGON ((-111.049728 44.488163, -111.050245 4... 29.6
Washington USA-3519 6724540 POLYGON ((-116.998073 46.33017, -116.906528 46... 27.2

现在我们可以绘制它:

import geoplot.crs as gcrs
gplt.cartogram(result, scale='Percent', projection=gcrs.AlbersEqualArea())
<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x11c3bda20>
../_images/Working_with_Geospatial_Data_36_1.png

保存格式

可以使用从地理空间文件格式中读取数据 GeoDataFrame.from_file . 可以使用将数据写入地理空间文件格式 GeoDataFrame.to_file . 默认情况下,这些方法将推断文件格式并保存到 Shapefile 分别是。要指定显式文件格式,请将该格式的名称传递给 driver 争论。例如:

nyc_boroughs.to_file('boroughs.geojson', driver='GeoJSON')

地理空间数据最简单且越来越常见的保存格式是 GeoJSON . geojson文件可能具有 .geojson.json 扩展,并以可读格式存储数据:

{
  "type": "Feature",
  "geometry": {
    "type": "Point",
    "coordinates": [125.6, 10.1]
  },
  "properties": {
    "name": "Dinagat Islands"
  }
}

从历史上讲,最常见的地理空间数据格式是 Shapefile . shapefile实际上不是文件,而是文件夹或文件夹中的文件组 zip 将这些数据一起存档可以对有关数据的非常复杂的信息进行编码。shapefile是一种二进制文件格式,因此它们不像GeoJSON文件那样是人类可读的,但是可以有效地对过于复杂的数据进行编码,以便于在GeoJSON中存储。

这是两种最著名的文件格式,但是 many many others . 获取支持的地理空间文件格式列表 geopandas refer to the fiona user manual .