1  GeoPandas介绍

        GeoPandas是一个开源项目，可以更轻松地使用python处理地理空间数据。

        GeoPandas扩展了Pandas中使用的数据类型DataFrame，允许对几何类型进行空间操作。

        GeoPandas的目标是使在python中使用地理空间数据更容易。它结合了Pandas和Shapely的能力，提供了Pandas的地理空间操作和多种Shapely的高级接口。GeoPandas可以让您轻松地在python中进行操作，否则将需要空间数据库，如PostGIS。

1.1 GeoPandas 概念

        GeoPandas 的核心是 geopandas.GeoDataFrame， pandas.DataFrame的一个子类。这个GeoDataFrame可以存储几何列并执行空间操作。

  几何图形由geopandas.GeoSeries处理。它是panda.series的一个子类。   

  GeoDataFrame 是一组Series的组合（可以是数字、文本、布尔变量等等），同时GeoDataFrame也包括了带有几何图形信息的geopandas.GeoSeries（点、多边形等）

   

         每个GeoSeries可以包含任何几何类型，并具有一个GeoSeries.crs属性，它存储关于投影的信息(crs代表坐标参考系统　Coordinate Reference System)。因此，GeoDataFrame中的每个GeoSeries可以在不同的投影中，例如，允许您拥有相同几何图形的多个版本(不同投影)。

        GeoDataFrame中只有一个GeoSeries被认为是活跃的几何图形(geometry所在的列），这意味着应用于GeoDataFrame的所有几何操作都是在这个活动列上操作的。

1.2 Geospatial数据的数据类型

• 地理空间数据使用一组不同的数据类型进行分析——> 使用 shapely 模块
  •  shapely 有一个名为 geometry 的类，其中包含不同的几何对象。 

 1.2.1 Point

点是表示二维空间中单个位置的对象( XY 坐标)。 

p1 = Point(0,0)
print(p1)
'''
POINT (0 0)
'''

• 注意，打印 p1 时，输出是 POINT (0 0)。 这表明返回的对象不 Python 的内置数据类型。

• print(p1 == (0,0))
  '''
  False
  '''
  
  
  print(type(p1))
  '''
  <class 'shapely.geometry.point.Point'>
  '''

1.2.1.2 Point转列表 

for i in states.head().geometry:
    print(list(i.coords))
'''
[(103.96085701087, 1.34617500824454, 0.0)]
[(103.978978976056, 1.38990401870983, 0.0)]
[(103.923130995416, 1.4099130157165, 0.0)]
[(103.866966992109, 1.37923700487091, 0.0)]
[(103.829073008381, 1.37276100286023, 0.0)]
'''

1.2.1.3 buffer 

圆

from shapely.geometry import Point
pt=Point(0,0)
patch=pt.buffer(10.0)
patch

 

 

1.2.2  Polygon

• 多边形，一个二维表面（一系列点）。
• 因为多边形由多个点组成，所以 shapely 多边形对象将元组列表作为参数。

• from shapely.geometry import Polygon
  polygon = Polygon([(0,0),(1,1),(1,0)])
  print(polygon)
  '''
  POLYGON ((0 0, 1 1, 1 0, 0 0))
  '''
  

1.2.3 LineString

一条直线

from shapely.geometry import LineString

line=LineString([(0,1),(2,2)])
line

 

２读取和写文件

 2.1 读文件

        假设我们有一个文件，这个文件同时有数据和几何信息（比如，GeoPackage, GeoJSON，Shapefile等），我们可以用geopandas.read_file来读取它。这个操作可以自动探测文件属性，并创建一个GeoDataFrame。

        比如我们使用’nybb‘数据集，这个数据集是纽约各区地图，GeoPandas安装的一部分。

import geopandas

path_to_data = geopandas.datasets.get_path("nybb")
gdf = geopandas.read_file(path_to_data)

gdf

2.1.1 GeoJSON 

• GeoJSON 是一种表示地理对象的格式。
• 它与常规 JSON 不同，因为它支持几何类型，例如：Point、LineString、Polygon、MultiPoint、MultiLineString、MultiPolygon 和 GeometryCollection。
• 使用 GeoJSON，使可视化变得更加容易， 这主要是因为 GeoJSON 允许存储几何数据类型

以Parks-Data.gov.sg为例：

import geopandas as gpd
states = gpd.read_file('parks/parks-geojson.geojson')
print(states.head())

'''
    Name                                        Description  \
0  kml_1  <center><table><tr><th colspan='2' align='cent...   
1  kml_2  <center><table><tr><th colspan='2' align='cent...   
2  kml_3  <center><table><tr><th colspan='2' align='cent...   
3  kml_4  <center><table><tr><th colspan='2' align='cent...   
4  kml_5  <center><table><tr><th colspan='2' align='cent...   

                              geometry  
0  POINT Z (103.96086 1.34618 0.00000)  
1  POINT Z (103.97898 1.38990 0.00000)  
2  POINT Z (103.92313 1.40991 0.00000)  
3  POINT Z (103.86697 1.37924 0.00000)  
4  POINT Z (103.82907 1.37276 0.00000)  
'''

2.1.2 KMX /KML

数据集以Singapore Police Force NPC Boundary-Data.gov.sg 为例

2.1.2.1 读取

import geotable
dir='C:/Users/16000/Downloads/singapore-police-force-npc-boundary/'
t = geotable.load(dir+'spf-boundaries.kml') #kml同样可以
t
'''
ame	NPC_NAME	DIVISION	DIV	INC_CRC	FMEL_UPD_D	geometry_object	geometry_layer	geometry_proj4
0	kml_1	Bukit Merah East Neighbourhood Police Centre	Central Divisional HQ	A	E2C8B705E260A602	20220113151514	POLYGON Z ((103.835389913422 1.29238449373637 ...	NPC_BOUNDARY	+proj=longlat +datum=WGS84 +no_defs
1	kml_2	Jurong East Neighbourhood Police Centre	Clementi Divisional HQ	D	D700372C36BA2227	20220113151514	POLYGON Z ((103.670422201851 1.23462116894326 ...	NPC_BOUNDARY	+proj=longlat +datum=WGS84 +no_defs
2	kml_3	Marina Bay Neighbourhood Police Centre	Central Divisional HQ	A	CF58A6E322230247	20220113151514	POLYGON Z ((103.867887287915 1.30391238416607 ...	NPC_BOUNDARY	+proj=longlat +datum=WGS84 +no_defs
3	kml_4	Jurong West Neighbourhood Police Centre	Jurong Divisional HQ	J	62E89DB5EB7BE0E9	20220113151514	POLYGON Z ((103.694482039648 1.30523552291172 ...	NPC_BOUNDARY	+proj=longlat +datum=WGS84 +no_defs
'''

2.2 写文件 

写文件直接用GeoDataFrame.to_file即可

默认的文件格式是Shapefile，但是你也可以用driver来赋值给

3 简单的方法 

        现在我们有了GeoDataFrame，可以开始处理它的几何属性了。

        因为我们只从文件中读取了一个几何列，它被自动视为活跃几何属性。GeoDataFrame上定义的方法将应用到“几何”列。

3.1 测量面积

        计算一个多边形的面积，我们可以使用GeoDataFrame.area属性，他会返回一个 pandas.Series。

   注意GeoDataFrame.area只会应用于活跃的几何列

gdf = gdf.set_index("BoroName")
gdf

gdf["area"] = gdf.area
gdf

 3.2 得到多边形的边界

 GeoDataFrame.boundary.

 3.3 得到多边形的质心

3.4 测量距离 

我们还可以测量每个质心与第一个质心的距离。

 这里我们先设定一个 Point作为我们的基准点，然后对某一列进行求距离的操作

        求得的结果是一个DataFrame，所以我们可以在地理空间数据集上使用所有的pandas功能，并使用属性和几何信息一起做数据操作。

3.5 制作地图

        GeoPandas还可以绘制地图，这样我们就可以检查我们的几何图形在空间中的样子。

        方法就是GeoDataFrame.plot()

        在下面的例子中，我们对gdf画图

3.5.1  首先我们查看当前的活跃几何列

       

3.5.2 然后绘图 

（这个相当于底图）

 我们用当前的活跃几何列描绘某一列GeoSeries

 3.5.3 切换活跃列

 我们切换活跃列（set_geometry），出来的就是不一样的东西

     

 3.5.4 叠加GeoSeries

我们也可以将两个GeoSeries叠加在一起。我们只需要用一个图作为另一个图的轴axis。

 

 gdf["centroid"]是在ax之上画成的，所以ax是“最底层”。

因而右边的图，ax交换之后，就显示不出黑点了，因为gdf["centroid"]被 gdf["geometry"]覆盖住了

ax=gpd.GeoSeries(t.geometries).plot(figsize=(25,7))
states['geometry'].plot(ax=ax,color='red')

3.6 Geopandas数据集

naturalearth_lowres
naturalearth_cities

 

 

4 创建几何图形

在我们已经有的几何图形的基础上，我们可以创建新的图形

4.1 凸包 convex hull

GeoDataFrame.convex_hull.

 4.2 缓冲 buffer

        在某些情况下，我们可能需要使用GeoDataFrame.buffer()缓冲几何图形。

        GeoDataFrame.buffer()方法会自动应用于活跃的几何列，但我们也可以直接应用于任何GeoSeries。

        让我们缓冲 区和它们的中心，并把它们画在一起。

缓冲活跃几何列10000英尺

gdf["buffered"] = gdf.buffer(10000)

 缓冲质心几何列10 000英尺

gdf["buffered_centroid"] = gdf["centroid"].buffer(10000)

做图可视化

ax = gdf["buffered"].plot(alpha=.5,figsize=(15,15))
gdf["buffered_centroid"].plot(ax=ax, color="red", alpha=.5)
gdf["bound"].plot(ax=ax, color="white", linewidth=.5)

 5 几何相关性

我们也可以去研究不同几何列的空间关系。利用上面的几何图形，我们可以检查哪些缓冲区与布鲁克林的原始几何体相交，也就是说，距离布鲁克林10000英尺以内的区域。

首先我们先找到布鲁克林区的多边形

gdf.loc["Brooklyn", "geometry"]

5.1 相交

然后我们可以看哪些 gdf["buffered"] 里面的几何体和布鲁克林区相交

 5.2 在原始多边形内

我们可以检查哪些缓冲质心是完全在原多边形内的。

 画图认证:

gdf = gdf.set_geometry("buffered_centroid")
ax = gdf.plot("within", legend=True, categorical=True, legend_kwds={'loc': "upper left"}) 
# using categorical plot and setting the position of the legend
gdf["bound"].plot(ax=ax, color="black", linewidth=.5)  
# passing the first plot and setting linewitdth to 0.5

6 投影 

        每一个GeoSeries 都有一个坐标参考系 Coordinate Reference System (CRS)，GeoSeries.crs

        CRS告诉GeoPandas几何图形的坐标在地球上的位置。

        在某些情况下，CRS是地理坐标，这意味着坐标以纬度和经度表示。在这些情况下，其CRS是WGS84，其授权代码是EPSG:4326。

        让我们看看纽约区GeoDataFrame的投影

 使用英尺坐标的几何图形的授权代码是 EPSG:2263。我们可以通过GeoSeries.to_crs()将GeoSeries投影到别的授权编码上 

gdf = gdf.set_geometry("geometry")
boroughs_4326 = gdf.to_crs("EPSG:4326")
boroughs_4326.crs

两个编码的区别在于，我们之前两个点的举例是 120 000 - 280 000 (feet)，现在是 40.5 - 40.9 (degrees) 

7 GeoSeries

类似于pandas中的Series

7.1 创建

import geopandas
from shapely.geometry import Polygon
s=geopandas.GeoSeries([
    Polygon([(0,0),(2,2),(0,2)]),
    Polygon([(0,0),(1,1),(0,1)])
])
s

 

7.2 主要函数

area	计算面积
intersection	求交集
plot	画图

 

参考资料：Introduction to GeoPandas — GeoPandas 0.9.0 documentation