Python+GIS ≈Geopandas?

最近读《红楼梦》，读到了宁荣二府受元妃之命前往清虚观打醮，那人马纷纷、轿车云云的场面，好不气派！在清虚观里张法官有意与宝玉说亲，但是心里只有黛玉的宝玉自然反感至极。至第二日，宝玉因张道士的说亲心中大不自在，且因为黛玉中暑，闷闷不乐。为此专门前往潇湘馆来看黛玉，却不知原为好意的探视却因为黛玉的偏执、任性的小性子惹起了两人的误会。至此大观园里的潇湘馆两人又是哭又是闹，甚至惊动了贾母。两人虽各处一室、却情发一心，一个掩面啜泣，一个仰天长叹，真真叫人看的捉急。

也难怪，哪个少男少女谈恋爱不是在这样的感情里互相较量过来的呢！

正所谓：

颦颦眉间泛情思，

冉冉身影惹人怜。

之前有很长一段时间觉得Python作为一种弱类型的语言，写起来很不得劲儿，写法总觉得怪怪的，可能是强类型的语言写惯了或者受Java亦或Golang先入为主的影响。不太想去接触Python，觉得语法又慢又不“正宗”，即便它很强大，解决问题很锋利，但依然让我觉得这不是一门值得学习的编程语言。

直到我看了一篇帖子，用几行代码就解决了一个问题，刹那间，好像突然就改变了我对Python的看法。

其实语言作为一种工具，并无优劣之分，不管是静态还是动态语言，都是站在编译器的肩膀上做一些解决业务问题的工作。Python作为高级语言的一种，简洁的语法、强大的第三方库这是它最大的优势，只要能给我们解决问题带来便利，它就是好工具。

今天我们来一起学习一下，强大的geopandas库。

Front

因为使用pip安装很多第三方库受限服务器都在国外的原因，所以下载速度很慢。

因此为了加速下载，我们可以配置全局的pip镜像下载源，步骤如下：

1. 使用Python的配置步骤：

  • 使用win+E打开文件管理器，在地址栏键入%appdata%；

  • 在该目录下新建pip文件夹；

  • 在pip文件夹下，新建pip.ini文件；

  • 在pip.ini文件内，输入以下内容：

[global]

timeout=6000 index-url = http://mirrors.aliyun.com/pypi/simple

trusted-host = mirrors.aliyun.com

至此，我们可以流畅的使用pip安装任何第三方库了。

2. 使用Anaconda的配置步骤：

  • 进入Anaconda安装目录的Scripts目录；

  • 右键启动cmd命令行窗口，输入以下命令，并回车；

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/ #设置搜索时显示通道地址 conda config --set show_channel_urls yes

3. 查看配置结果；

conda config --show channels

[图片上传失败...(image-3cf46e-1648374378608)]

About GeoPandas

GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and matplotlib for plotting.

geopandas包详情:

[图片上传失败...(image-a61ecc-1648374378608)]

顾名思义，GeoPandas 通过添加对地理空间数据的支持来扩展流行的数据科学库 pandas。

GeoPandas 中的核心数据结构是 geopandas.GeoDataFrame，它是pandas.DataFrame 的子类，可以存储几何列并执行空间操作。 geopandas.GeoSeries 是 pandas.Series 的子类，用于处理几何图形。因此，您的 GeoDataFrame是pandas.Series 与传统数据（数字、布尔值、文本等）和 geopandas.GeoSeries 与几何（点、多边形等）的组合。

Install

如果你使用Anaconda,那么键入以下命令即可：

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">conda install geopandas </pre>

如果你是使用Python，那么首先去下面的网站下载对应的依赖，因为pip不会帮助我们自己安装依赖。

https://www.lfd.uci.edu/~gohlke/pythonlibs/

依赖包根据自己的python版本和OS选择合适的即可，包括

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">Fiona pyproj rtrre shapely GDAL </pre>

实际按照上面的顺序安装会报错，因此实际的安装顺序如下：

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">GDAL->Fiona->pyproj->shapely->rtree </pre>

安装命令为,以GDAL为例，其它亦然：

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">pip install GDAL-3.4.2-cp37-cp37m-win_amd64.whl </pre>

验证一下是否安装成功：

[图片上传失败...(image-cc151e-1648374378607)]

Example

注：以下代码均可单独运行，并输出结果至控制台

geopandas.GeoSeries

构造POINT

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from shapely.geometry import Point
from geopandas import geopandas

s = geopandas.GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)])
print(s)` </pre>

构造具有坐标系的POINT

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from shapely.geometry import Point
from geopandas import geopandas

s = geopandas.GeoSeries(
[Point(1, 1), Point(2, 2), Point(3, 3)], crs="EPSG:3857")` </pre>

geopandas.GeoSeries.area

1. 构造地理集合，并量算面积

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">from shapely.geometry import Polygon, LineString, Point s = geopandas.GeoSeries( [ Polygon([(0, 0), (1, 1), (0, 1)]), Polygon([(10, 0), (10, 5), (0, 0)]), Polygon([(0, 0), (2, 2), (2, 0)]), LineString([(0, 0), (1, 1), (0, 1)]), Point(0, 1) ] ) area = s.area print(area) </pre>

geopandas.GeoSeries.boundary

构造地理集合，并返回低维对象(抽稀)

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from shapely.geometry import Polygon, LineString, Point

from geopandas import geopandas
s = geopandas.GeoSeries(
[
Polygon([(108.923166,35.363084), (109.322899,35.274691), (109.355755,35.283365),(108.936301,35.380966),(108.923166,35.363084)]),
LineString([(0, 0), (1, 1), (1, 0)]),
Point(0, 0),
]
)
boundary = s.boundary
print(boundary)` </pre>

geopandas.GeoSeries.x

获取地理对象的X值

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from shapely.geometry import Point
from geopandas import geopandas

s = geopandas.GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)])
s.x
print(s.x)` </pre>

geopandas.GeoSeries.within

返回两个几何图形是否包含，值为bool

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from shapely.geometry import Polygon, LineString, Point
from geopandas import geopandas

s1 = geopandas.GeoSeries(
[
Polygon([(0, 0), (2, 2), (0, 2)]),
Polygon([(0, 0), (1, 2), (0, 2)]),
LineString([(0, 0), (0, 2)]),
Point(0, 1),
],
)
s2 = geopandas.GeoSeries(
[
Polygon([(0, 0), (1, 1), (0, 1)]),
LineString([(0, 0), (0, 2)]),
LineString([(0, 0), (0, 1)]),
Point(0, 1),
],
index=range(1, 5),
)

polygon = Polygon([(0, 0), (2, 2), (0, 2)])
within = s1.within(polygon)

print(within)

print(s1.within(s2))

print(s1.within(polygon))` </pre>

geopandas.GeoSeries.to_json

地理集合转json

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">from shapely.geometry import Point s = geopandas.GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)]) s.to_json() </pre>

geopandas.GeoDataFrame.set_crs

GeoDataFrame设置CRS

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`import os
import sys
import pyproj
from shapely.geometry import Point
from geopandas import geopandas

os.environ['PROJ_LIB'] = os.path.dirname(sys.argv[0])

os.environ['PROJ_LIB'] = "D:\ProgramData\Anaconda3\Library\share\proj\"

d = {'col1': ['name1', 'name2'], 'geometry': [Point(1, 2), Point(2, 1)]}
gdf = geopandas.GeoDataFrame(d)
print(gdf.crs)
gdf = gdf.set_crs('epsg:4490')
print(gdf.crs)` </pre>

geopandas.read_file

读取矢量文件

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from geopandas import geopandas
from geopandas.geodataframe import GeoDataFrame

df = geopandas.GeoDataFrame.from_file("W:\desktop\葛瑶\毕业论文\GIS\shp\城六区遗址公园面状.shp")

head = df.head(6)
print(head)` </pre>

Missing and empty geometries

缺省和空值

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`from shapely.geometry import Polygon

from geopandas import geopandas

s = geopandas.GeoSeries([Polygon([(0, 0), (1, 1), (0, 1)]), None, Polygon([])])

print(s.area)
print(s.union(Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])))
print(s.intersection(Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])))` </pre>

1. Creating a GeoDataFrame from a DataFrame with coordinates

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`import pandas as pd
import geopandas
import matplotlib.pyplot as plt

df = pd.DataFrame(
{'City': ['Buenos Aires', 'Brasilia', 'Santiago', 'Bogota', 'Caracas'],
'Country': ['Argentina', 'Brazil', 'Chile', 'Colombia', 'Venezuela'],
'Latitude': [-34.58, -15.78, -33.45, 4.60, 10.48],
'Longitude': [-58.66, -47.91, -70.66, -74.08, -66.86]})

gdf = geopandas.GeoDataFrame(
df, geometry=geopandas.points_from_xy(df.Longitude, df.Latitude))

print(gdf.head())

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

We restrict to South America.

ax = world[world.continent == 'South America'].plot(
color='white', edgecolor='black')

We can now plot our GeoDataFrame.

gdf.plot(ax=ax, color='red')

plt.show()` </pre>

运行结果如下：

[图片上传失败...(image-6093b0-1648374378607)]

2. create a Isodensity map

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">`import pandas as pd
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point, geo # 经纬度转换为点
import adjustText as aT
import mapclassify
import warnings

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

map_name = "西安市区域专题图"

warnings.filterwarnings('ignore')

读取数据为geodataframe格式

geo_ = gpd.GeoDataFrame.from_file('W:\Google_Download\西安市.json', encoding='utf-8')
print(geo_)
geo_.head(3)
df = pd.read_csv(r'C:\Users\王斌\Desktop\POI.csv', encoding='utf-8')

df['geometry'] = list(zip(df['wgs84_lon'], df['wgs84_lat'])) # 经纬度组合为新列geometry，与形状里的该列对应
df['geometry'] = df['geometry'].apply(Point) # 经纬度转化为坐标点
gpd_df = gpd.GeoDataFrame(df)

base = geo_.plot(color='lightblue', edgecolor='grey', figsize=(15, 15)) # 利用形状信息画底图
gpd_df.plot(ax=base, color='red', marker="o", markersize=50, alpha=0.01) # 在底图上添加出租房位置
plt.gca().xaxis.set_major_locator(plt.NullLocator()) # 去掉x轴刻度
plt.gca().yaxis.set_major_locator(plt.NullLocator()) # 去掉y轴刻度

geo_.plot(figsize=(8, 6), # 图像大小
column='name', # 分级设色字段
# scheme='quantiles', # MapClassify-分级类型
legend=True, # 图例
cmap='Pastel1_r', # 渐变色带的名称#Set2
edgecolor='k') # 边框颜色

图名

plt.Artist

plt.title('{}'.format(map_name), fontsize=18, fontweight='bold')

显示网格，透明度为50%

plt.grid(True, alpha=0.5)

plt.show()

plt.savefig('./image/{}.png'.format(map_name), dpi=300)` </pre>

如下为运行结果，红色斑块为西安市城六区的POI点，为21万余点。

[图片上传失败...(image-be99a2-1648374378607)]

[图片上传失败...(image-c95103-1648374378607)]

Summary

geopandas的函数较多，这里我选了几个运行和学习，大家有兴趣可以根据自己的需要进行学习。

geopandas作为一个Python与GIS结合的库，根据官网的介绍其提供的功能还是比较丰富的，能够进行较多的地理数据处理，也可以快速绘制一些基础地图，并能够灵活的对地图进行设计。

我感觉可以满足大部分的GIS小场景，但是遇到比较复杂的业务场景或者功能，可能还得借助桌面端的软件来辅助完成。

总体来讲，我感觉在解决一些小的GIS问题上来说还是很锋利的。

值得一学！

Error Summary

这里额外记录一个错误，我排查、上网找资料好久才算解决。

这里记录下来供大家参考！

错误的原因就是使用geopandas读取矢量文件或者绘图时会引发pyproj库的无效投影错误，错误截图如下：

pyproj.ecception.CRSError: Invalid projection:epsg:2263:(Internal Proj Error:···)

[图片上传失败...(image-7da2cd-1648374378607)]

这是在Proj官网上给出的一个常见错误排查方案，能提供思路，但是没有解决问题。

Why am I getting the error “Cannot find proj.db”?

The file proj.db must be readable for the library to properly function. Like other resource files, it is located using a set of search paths. In most cases, the following paths are checked in order:

  • A path provided by the environment variable PROJ_LIB.
  • A path built into PROJ as its resource installation directory (typically ../share/proj relative to the PROJ library).
  • The current directory.

Note that if you’re using conda, activating an environment sets PROJ_LIB to a resource directory located in that environment.

最后在根据在网上查看许久发现线索指向Proj.db，该文件可能有问题，因为报错总是提示SQLite error,具体来说是关于初始化查询的一个错误，说明 Proj.db可能是有问题的。

使用如下代码打印出目前pyproj的资源指向路径：

<pre class="custom" data-tool="mdnice编辑器" style="margin-top: 10px; margin-bottom: 10px; border-radius: 5px; box-shadow: rgba(0, 0, 0, 0.55) 0px 2px 10px;">data_dir = pyproj.datadir.get_data_dir() print(data_dir) </pre>

在我的笔记本上该路径指向D:\ProgramData\Anaconda3\Lib\site-packages\pyproj\proj_dir\share\proj

经查看该文件下的Proj.db为7.7M左右，但是pyproj官方提示的在Anaconda的安装目录/share/proj下的Proj.db大小为5.36M。

实在想不到解决办法的我，大胆地把D:\ProgramData\Anaconda3\Lib\site-packages\pyproj\proj_dir\share\proj下的Proj.db替换为了Anaconda的安装目录/share/proj下的Proj.db。

正所谓：“无心插柳柳成荫”，没想到问题好像解决了。程序可以正确的进行矢量文件读取和制图。

按道理来讲，由Anaconda下载的依赖，其中的依赖关系会被较好的处理，很难出现这样的错误。

我记得我在网上好像看到了Pyproj库是建立在EPSG数据库之上的，这个错误可能是由与版本冲突造成的吧。

希望后面可以搞清楚原因。

知道其中缘由的同学，一定要私信我哈~

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,362评论 5赞 477
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,330评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,247评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,560评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,580评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,569评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,929评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,587评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,840评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,596评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,678评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,366评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,945评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,929评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,165评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 43,271评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,403评论 2赞 342

推荐阅读更多精彩内容