2023年3月

本文主要介绍GeoPandas的基本使用方法,以绘制简单的地图。GeoPandas是一个Python开源项目,旨在提供丰富而简单的地理空间数据处理接口。GeoPandas扩展了Pandas的数据类型,并使用matplotlib进行绘图。GeoPandas官方仓库地址为:
GeoPandas
。GeoPandas的官方文档地址为:
GeoPandas-doc

本文所有代码见:
Python-Study-Notes

GeoPandas推荐使用Python3.7版本及以上,运行环境最好是linux系统。GeoPandas安装命令如下:

pip install geopandas

如果上述命令安装出问题,则推荐使用conda安装GeoPandas,命令如下:

conda install geopandas

或:

conda install --channel conda-forge geopandas

# jupyter notebook环境去除warning
import warnings
warnings.filterwarnings("ignore")

# 查看geopandas版本
import geopandas as gpd

gpd.__version__

'0.10.2'

1 基础结构

GeoPandas在Pandas的基础上增加对地理空间数据的支持。其主要结构为GeoSeries和GeoDataFrame。如下图所示,整张图片表示一个GeoDataFrame,行表示地图上的某一区域,列表示该区域的属性数据,其中geometry表示该区域的位置信息。每一个区域的位置形状用一个或多个几何对象表示,几何对象在GeoPandas中以GeoSeries格式数据构成。

1.1 GeoSeries

1.1.1 GeoSeries基础元素

在Geogeopandas中,GeoSeries用于构成基础的几何对象。GeoSeries由Python库Shapely实现的几何对象构成,这些几何元素可以是点(集)、线(集)、面(集):

  • Points / Multi-Points

  • Lines / Multi-Lines

  • Polygons / Multi-Polygons

在jupyter notebook中, 可以直接以svg格式展示GeoSeries中的单个元素,使用如下:

Points 点

# Point(x,y)用于创建单个点
from shapely.geometry import Point
p = gpd.GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)],index=['p1','p2','p3'])
p
p1    POINT (1.00000 1.00000)
p2    POINT (2.00000 2.00000)
p3    POINT (3.00000 3.00000)
dtype: geometry
# 展示单个元素
p[0]

svg

Multi-Points 点集

# MultiPoint([(x1,y1), ... ,(xn,yn)])用于创建多个点构成的集合
from shapely.geometry import MultiPoint
mp = gpd.GeoSeries([MultiPoint([(0, 1), (1, 0)]),
                    MultiPoint([(0, 0), (1, 1)])],
                    index=['mp1', 'mp2'])
mp
mp1    MULTIPOINT (0.00000 1.00000, 1.00000 0.00000)
mp2    MULTIPOINT (0.00000 0.00000, 1.00000 1.00000)
dtype: geometry
# 展示单个元素
mp[0]

svg

Lines 线

# LineString([(x1,y1), ... , (xn,yn)])用于创建单条线段
from shapely.geometry import LineString
l = gpd.GeoSeries([LineString([(0.5, 0), (1, 2), (0, 1)]),
                   LineString([(0, 0), (3.3, 1)])],
                   index=['l1', 'l2'])
l
l1    LINESTRING (0.50000 0.00000, 1.00000 2.00000, ...
l2        LINESTRING (0.00000 0.00000, 3.30000 1.00000)
dtype: geometry
# 展示单个元素
l[0]

svg

MultiLineString 线集

# MultiLineString([LineString1,LineString2])用于创建多条线段的集合
from shapely.geometry import MultiLineString
ml = gpd.GeoSeries([MultiLineString([[(1, 0), (2, 1)], [(-1, 0), (0, 5), (1, 0)]]),
               MultiLineString([[(0.5, 0), (1, 2), (0, 1)]])],
               index=['m1','m2'])
ml
m1    MULTILINESTRING ((1.00000 0.00000, 2.00000 1.0...
m2    MULTILINESTRING ((0.50000 0.00000, 1.00000 2.0...
dtype: geometry
# 展示单个元素
ml[0]

svg

Polygon 面

在Polygon元素,其创建函数为:

# shell表示单个多边形坐标,hole表示是否在内部创建空洞
Polygon(shell=None, holes=None)
# 创建单个多边形
from shapely.geometry import Polygon

pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)]),
                     Polygon([(0, 0), (1, 1), (1, 2)])],
                    index=['pol1','pol2'])
pol
pol1    POLYGON ((0.00000 0.00000, 4.00000 1.00000, 5....
pol2    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 1....
dtype: geometry
# 展示单个元素
pol[0]

svg

# 创建带空洞的单个多边形
pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)],
                     [((1, 1),(1,2),(2,2),(2,1)),
                     ((4,2),(4,3),(1,5))])])
pol
0    POLYGON ((0.00000 0.00000, 4.00000 1.00000, 5....
dtype: geometry
# 展示单个元素
pol[0]

svg

MultiPolygon 面集

# MultiPolygon([Polygon,Polygon])用于创建多个多边形的集合
from shapely.geometry import Polygon, MultiPolygon

mpol = gpd.GeoSeries([MultiPolygon([Polygon([(0, 0), (0, 1), (1, 1), (1, 0)]),
                                    Polygon([(2, 3), (3, 2), (3, 3)])])])
mpol
0    MULTIPOLYGON (((0.00000 0.00000, 0.00000 1.000...
dtype: geometry
# 展示单个元素
mpol[0]

svg

1.1.2 GeoSeries属性与方法

GeoSeries常用属性如下:

  • area: 计算GeoSeries每一个元素的面积,面积单位取决于坐标系
  • bounds: 获得GeoSeries每一个元素的最小外接矩形的左上角坐标和右下角坐标
  • total_bounds:计算整个GeoSeries对象的最小外接矩形的左上角坐标和右下角坐标
  • length:计算GeoSeries每一个元素的周长,周长单位取决于坐标系
  • geom_type: 获得GeoSeries每一个元素的类型
  • is_valid: 判断GeoSeries每一个元素是否构成合理的多边形
  • centroid:获得GeoSeries每一个元素的几何中心
# 创建包含两个多边形的GeoSeries
from shapely.geometry import Polygon

pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)]),
                    Polygon([(0, 0), (1, 1), (1, 8)])],
                    index=['pol1','pol2'])
# 计算各个元素的面积
pol.area
pol1    21.0
pol2     3.5
dtype: float64
# 计算各元素的最小外接矩形的顶点数据
pol.bounds
minx miny maxx maxy
pol1 0.0 0.0 5.0 7.0
pol2 0.0 0.0 1.0 8.0
# 计算整个对象的最小外接矩形的顶点数据
pol.total_bounds
array([0., 0., 5., 8.])
# 计算GeoSeries每一个元素的周长
pol.length
pol1    19.762298
pol2    16.476471
dtype: float64
# 获得GeoSeries每一个元素的类型
pol.type
pol1    Polygon
pol2    Polygon
dtype: object
# 判断GeoSeries每一个元素是否构成合理的多边形
pol.is_valid
pol1    True
pol2    True
dtype: bool
# 获得GeoSeries每一个元素的几何中心
pol.centroid
pol1    POINT (1.88889 3.00000)
pol2    POINT (0.66667 3.00000)
dtype: geometry

GeoSeries常用方法如下:

  • distance(): 计算GeoSeries每一个元素到某一几何体的距离,具体使用查看官方文档
  • representative_point():返回位于GeoSeries每一个元素中的坐标点,该点不会是几何中心。
  • to_crs():改变GeoSeries的坐标参考系
  • plot():绘制整个GeoSeries
# 创建包含两个多边形的GeoSeries
from shapely.geometry import Polygon

pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)]),
                    Polygon([(0, 0), (1, 1), (1, 8)])],
                    index=['pol1','pol2'])
# 计算距离
point = Point(1, 0)
pol.distance(point)
pol1    0.242536
pol2    0.707107
dtype: float64
# 返回位于GeoSeries每一个元素中的坐标点
pol.representative_point()
pol1    POINT (1.25000 5.00000)
pol2    POINT (0.78125 4.50000)
dtype: geometry
# 绘制整个GeoSeries
pol.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fbe00a24090>

1.2 GeoDataFrame

GeoDataFrame是包含GeoSeries的表格数据结构,GeoDataFrame最重要的属性是它始终具有一个表示地理状态信息的GeoSeries列,即geometry列。调用GeoDataFrame的方法或属性时,将始终作用于geometry列。

下面用一个世界地图的实例介绍GeoDataFrame常用函数的使用方法,具体函数的接口请查阅官方文档。

⚠⚠⚠注意该世界地图数据集来自geopandas自带数据集,部分数据数据可能缺失或错误,且数据集仅可用于学习不做其他用途

import geopandas as gpd
# 加载自带世界地图数据集
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# 查看数据集
# 人口,大洲,国名,国家缩写,ISO国家代码,gdp,地理位置数据
world.head()
pop_est continent name iso_a3 gdp_md_est geometry
0 920938 Oceania Fiji FJI 8374.0 MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1 53950935 Africa Tanzania TZA 150600.0 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2 603253 Africa W. Sahara ESH 906.5 POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3 35623680 North America Canada CAN 1674000.0 MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4 326625791 North America United States of America USA 18560000.0 MULTIPOLYGON (((-122.84000 49.00000, -120.0000...
# 绘制GeoDataFrame
# column选择为地图着色的数据列
# kind设置图片类型,某些类型会导致column参数不可用
# cmap设置matplotlib调色板
# legend设置是否显示图例
# figsize设置图片大小
# legend_kwds设置matplotlib的legend参数,该参数需要根据实际图例形式设置
world.plot(column="continent",kind="geo", 
            cmap="coolwarm",legend=True,figsize=(9,6), 
            legend_kwds={'loc': (0,1.02),'ncol':4})
<matplotlib.axes._subplots.AxesSubplot at 0x7fbe00551790>

png

# 获得设置地理位置信息的列名
world.geometry.name
# 修改地理位置信息的列名,一般不要修改
# world = world.rename(columns={'geometry': 'borders'}).set_geometry('borders')
# world.geometry.name
'geometry'
# 计算各个区域的几何中心,并新增列centroid_column
world['centroid_column'] = world.centroid
# 以centroid_column表示地理位置信息
world = world.set_geometry('centroid_column')
world.head()
pop_est continent name iso_a3 gdp_md_est geometry centroid_column
0 920938 Oceania Fiji FJI 8374.0 MULTIPOLYGON (((180.00000 -16.06713, 180.00000... POINT (163.85316 -17.31631)
1 53950935 Africa Tanzania TZA 150600.0 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982... POINT (34.75299 -6.25773)
2 603253 Africa W. Sahara ESH 906.5 POLYGON ((-8.66559 27.65643, -8.66512 27.58948... POINT (-12.13783 24.29117)
3 35623680 North America Canada CAN 1674000.0 MULTIPOLYGON (((-122.84000 49.00000, -122.9742... POINT (-98.14238 61.46908)
4 326625791 North America United States of America USA 18560000.0 MULTIPOLYGON (((-122.84000 49.00000, -120.0000... POINT (-112.59944 45.70563)
# 根据gdp绘图
world.plot(column="gdp_md_est",legend=True)
<matplotlib.axes._subplots.AxesSubplot at 0x7fbe00489190>

png

# 叠加几何中心绘图
ax = world["geometry"].plot()
world["centroid_column"].plot(ax=ax, color="red")
<matplotlib.axes._subplots.AxesSubplot at 0x7fbe002cfe90>

png

# 设置数据小数精度
gpd.options.display_precision = 1
world.head()
pop_est continent name iso_a3 gdp_md_est geometry centroid_column
0 920938 Oceania Fiji FJI 8374.0 MULTIPOLYGON (((180.0 -16.1, 180.0 -16.6, 179.... POINT (163.9 -17.3)
1 53950935 Africa Tanzania TZA 150600.0 POLYGON ((33.9 -1.0, 34.1 -1.1, 37.7 -3.1, 37.... POINT (34.8 -6.3)
2 603253 Africa W. Sahara ESH 906.5 POLYGON ((-8.7 27.7, -8.7 27.6, -8.7 27.4, -8.... POINT (-12.1 24.3)
3 35623680 North America Canada CAN 1674000.0 MULTIPOLYGON (((-122.8 49.0, -123.0 49.0, -124... POINT (-98.1 61.5)
4 326625791 North America United States of America USA 18560000.0 MULTIPOLYGON (((-122.8 49.0, -120.0 49.0, -117... POINT (-112.6 45.7)

此外,pandas的数据操作也可以用于geopandas。

# 查看数据
world.iloc[:4,[3,4]]
iso_a3 gdp_md_est
0 FJI 8374.0
1 TZA 150600.0
2 ESH 906.5
3 CAN 1674000.0
# 查看人口小于10万的国家地区
world.loc[world['pop_est']<1e5,['pop_est','name','gdp_md_est','geometry']]
pop_est name gdp_md_est geometry
20 2931 Falkland Is. 281.8 POLYGON ((-61.2 -51.9, -60.0 -51.2, -59.1 -51....
22 57713 Greenland 2173.0 POLYGON ((-46.8 82.6, -43.4 83.2, -39.9 83.2, ...
23 140 Fr. S. Antarctic Lands 16.0 POLYGON ((68.9 -48.6, 69.6 -48.9, 70.5 -49.1, ...
159 4050 Antarctica 810.0 MULTIPOLYGON (((-48.7 -78.0, -48.2 -78.0, -46....

2 数据集基本操作

2.1 数据集读写

对于包含地理信息数据的文件(例如GeoPackage、GeoJSON、Shapefile),geopandas提供了read_file()函数接口来自动读取对应的类型文件并返回GeoDataFrame对象。此外,geopandas也提供to_file()函数接口以实现不同类型的地理信息文件写入。下面以一个实际例子来介绍数据集的读写。

对于地理信息数据,通常的获取方式是从专业网站下载,或从谷歌地图和高德地图获取。计算机领域常用的地理信息数据文件格式为GeoJSON,关于GeoJSON格式的中国地理信息数据,推荐从阿里云数据可视化平台获取:
DataV.GeoAtlas
。DataV.GeoAtlas提供了我国各省各市区县地图信息(不包含乡镇)。

如下图所示,打开网页后,点击或查找所需的区域->选择需要的数据格式->选择是否包含子区域->复制JSON API->read_file()函数读取文件即可。例如本文选择江苏省的地理信息,包含子区域的意思为是否包含江苏省内地区的地理数据。

import geopandas as gpd

# 读取数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")
# 展示数据
# abcode:行政地区编码;
# childrenNum,子区域个数,在这里表示城市所辖区县数,
# parent:父区域编码,320000表示江苏省的意思
data.head()
adcode name childrenNum level parent subFeatureIndex geometry
0 320100 南京市 11 city {'adcode': 320000} 0 MULTIPOLYGON (((119.1 32.5, 119.1 32.5, 119.1 ...
1 320200 无锡市 7 city {'adcode': 320000} 1 MULTIPOLYGON (((119.5 31.2, 119.5 31.1, 119.6 ...
2 320300 徐州市 10 city {'adcode': 320000} 2 MULTIPOLYGON (((118.4 34.4, 118.4 34.4, 118.4 ...
3 320400 常州市 6 city {'adcode': 320000} 3 MULTIPOLYGON (((120.0 32.0, 120.0 32.0, 120.0 ...
4 320500 苏州市 9 city {'adcode': 320000} 4 MULTIPOLYGON (((119.9 31.2, 119.9 31.2, 119.9 ...
# 展示江苏省地图,横轴为经度,纵轴为纬度
data.plot(figsize=(6,6),color='green')
<matplotlib.axes._subplots.AxesSubplot at 0x7fbe002057d0>

png

# 写入Shapefile
data.to_file("countries.shp")

# 写入GeoJSON
# 默认格式是Shapefile,通过driver参数可以设置其他类型
data.to_file("countries.geojson", driver='GeoJSON')

2.2 数据集筛选

GeoPandas提供了各种数据集筛选函数,以提取感兴趣区域。其在0.1.0版本中新增了Bounding Box Filter,在0.7.0版本中新增了Geometry Filter和Row Filter,以及其他需要Fiona库的过滤器。

2.2.1 Bounding Box Filter

GeoPandas在0.1.0版本引入了Bounding Box Filter。Bounding Box Filter用于提取与边界框相交的子区域,注意边界框格式为(左下角x, 左下角y, 右上角x, 右上角y)。

# 框选数据
import geopandas as gpd

# 设置框选范围
bbox = (119, 32, 121, 33)
# 读取数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json", bbox=bbox)

# 结果是筛选与bbox相交的子区域
data.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fbe000e1c90>

png

# 查看边界框与哪些子区域相交
# 在地图上绘制边界框
from shapely import geometry

# 设置框选范围
bbox = (119, 32, 121, 33)
# 读取数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")
ax = data.plot()
# 在地图上绘制bbox框
ax = gpd.GeoSeries([geometry.box(minx=bbox[0],miny=bbox[1], 
                                 maxx=bbox[2],maxy=bbox[3]).boundary]).plot(ax=ax, color='red')

png

2.2.2 Geometry Filter

GeoPandas在0.7.0版本引入了Geometry Filter。与Bounding Box Filter类似,Geometry Filter也是筛选与指定区域相交的子区域,但Geometry Filter输入的是几何图形。

# 查看几何图形与哪些子区域相交
from shapely import geometry

# 设置框选范围
p = geometry.Polygon([(119,32),(120,33), (121, 33)])
# 读取数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json",mask=p)
# 筛选区域
data.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6fbf590>

png

# 查看几何图形与哪些子区域相交
from shapely import geometry

# 设置框选范围
p = geometry.Polygon([(119,32),(120,33), (121, 33)])
# 读取数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")

ax = data.plot()
ax = gpd.GeoSeries([p.boundary]).plot(ax=ax, color='red')

png

2.2.3 Row Filter

GeoPandas在0.7.0版本引入了Row Filter。Row Filter用于提取指定行范围的数据。

# 读取前3行数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json",rows=3)
data
adcode name childrenNum level parent subFeatureIndex geometry
0 320100 南京市 11 city {'adcode': 320000} 0 MULTIPOLYGON (((119.1 32.5, 119.1 32.5, 119.1 ...
1 320200 无锡市 7 city {'adcode': 320000} 1 MULTIPOLYGON (((119.5 31.2, 119.5 31.1, 119.6 ...
2 320300 徐州市 10 city {'adcode': 320000} 2 MULTIPOLYGON (((118.4 34.4, 118.4 34.4, 118.4 ...
# 读取第1行和第2行数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json",rows=slice(1, 3))
data
adcode name childrenNum level parent subFeatureIndex geometry
0 320200 无锡市 7 city {'adcode': 320000} 1 MULTIPOLYGON (((119.5 31.2, 119.5 31.1, 119.6 ...
1 320300 徐州市 10 city {'adcode': 320000} 2 MULTIPOLYGON (((118.4 34.4, 118.4 34.4, 118.4 ...

2.2.4 cx索引器

除了标准的Pandas方法外,GeoPandas还提供使用cx索引器进行基于坐标的索引。通过cx索引器可以挑选特定区域。

# 读取世界地图
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# 筛选与经度范围为80到200,纬度小于0相交的区域,这里坐标代表经纬度。
southern_world = world.cx[80:200, :0]

southern_world.plot(figsize=(9, 6))
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6f03150>

png

# 查看与哪些子区域相交
# 读取世界地图
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

ax = world.plot(color='grey',alpha=0.5)

# 筛选与经度范围为80到200,纬度小于0相交的区域,这里坐标代表经纬度。
southern_world = world.cx[80:200, :0]

southern_world.plot(ax=ax,figsize=(9, 6))
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6ee2450>

png

3 坐标参考系

3.1 坐标参考系简介

地信、计算机和遥感等领域的从业人员或多或少都会接触地理信息系统(GIS,Geographic Information System)的相关知识。所谓GIS简单来说就是一个以计算机为核心,对地理空间位置相关数据进行创建、管理、分析、绘制和展示的多功能集成信息系统。绘制地图,则需要了解GIS中的坐标参考系。本文只是简单介绍坐标参考系的相关内容。相关内容总结于以下文章,想要具体了解更多内容也可以看看这些文章。

如下图所示,对于二维平面上的任何一点,我们可以通过如下的笛卡尔坐标系确定其唯一坐标,并用
\((x,y)\)
表示。

然而我们生活在一个恰好是“圆形”的三维地球上,一般来说要定义地球上物体的位置,需要一个适应地球形状的三维坐标系。但当我们在纸上或计算机屏幕上展示地球时,受限于存储介质(纸张、沙盘、书本、屏幕)原因,不得不使用二维平面来表达三维的地球。这种将三维球面转换为二维平面所用到的坐标系参考就是坐标参考系统(CRS,Coordinate Reference System),如下图所示。

常见情况下,根据不同的转换方式,CRS又可分为地理坐标系和投影坐标系。

地理坐标系统

参考平面为球面,用球面坐标来表示地球上的位置,叫做地理坐标系统。这些坐标系的坐标单位通常为角度,如经纬度。所谓经度就是指离被称为本初子午线的南北方向走线以东或以西的度数。纬度则是指地球某点与球心的连线和地球赤道面所成的线面角。如下图所示:

在数据表示中,东经正数,西经为负数,经度范围为-180度到180度。赤道以北为北纬,以南为南纬,纬度范围为-90度到90度。按照常用地理坐标系统表示的二维地图如下所示:

投影坐标系统

由于我们所在的地球实际不是规则球体,一经度和一纬度代表的实际距离是不相等的(赤道除外)。这导致我们无法通过经纬度坐标精确计算某块地区的面积,长宽。也无法精准表达某一块地区的实际形状。因此,有必要采取某些方法将地球表面投影到一个二维平面,这种地理投影方法也就是投影坐标系,实际作用为把地球表面经纬度按照一定的数学规则转换为平面坐标。

投影坐标系在二维平面中定义,其始终基于地理坐标系,位置由二维坐标标识,度量单位用米、千米等来表示。不同的投影方法会导致不同的投影结果,因此各种各样的投影坐标系统被提出。但是不管是哪种地图投影方法,都会导致地图投影变形。三维投影二维出现投影变形是不可避免的,但不会在任何位置都会产生变形。那些没有发生变形的点或线称为标准点线,距离标准点线越近,变形程度则越小。因此在实际投影某块区域时,最好选择标标准点线离该区域近的投影坐标系。

关于常见投影方法的介绍可以阅读本节开头给出的参考文章。想要感受不同投影方法带来的效果,可以使用一个非常酷的坐标系投影介绍网站:
worldmapcreator
。在worldmapcreator中,我们可以选择不同的投影方法来实时显示投影变换效果,以及自定义投影结果,具体如下:

3.2 EPSG编码

EPSG:European Petroleum Survey Group(欧洲石油调查小组)是一个涉及测地学、测量、制图学与石油勘探相关的科学组织,它维护和发布的EPSG编码系统提供了各种坐标参照系统crs的数据集参数。我们通过如下网站查询EPSG投影坐标的编码:
espg.io
,该网站界面介绍如下:

互联网常见的坐标参考系如下图所示。WGS84和CGCS2000都为地理坐标系统。其中WGS84(World Geodetic System 1984)是目前应用最为广泛的坐标系统。
在没有特别声明的时候,默认使用WGS84坐标系统存储地图数据或提供地图数据
。WGS84坐标单位为度,WGS84的EPSG代码为EPSG:4326。CGCS2000是我国当前最新的国家大地坐标系,它的EPSG编码为EPSG:4490,在大多数精度要求不高的情况下,我们可以默认WGS84坐标系等同于CGCS2000坐标系。

除此之外,WGS84 Web Mercator(EPSG:3857)坐标系统也较为常用。EPSG:3857为投影坐标系统,单位为米。EPSG:3857在WGS84坐标系基础上进行伪墨卡托投影(Pseudo-Mercator)得到。关于Pseudo-Mercator介绍见
墨卡托投影
。墨卡托投影能够保证方向的正确,所以EPSG:3857广泛用于导航海航。但越到高纬度,墨卡托投影会导致大小扭曲越严重,到两极会被放到无限大。在长度表示上EPSG:3857有很大偏差,无法用于实际距离和面积的测算。墨卡托投影下每个国家的大小和实际大小的差异见下图:

当然CRS还有其他常见的数据编码格式。计算机领域更多地使用EPSG编码。

3.3 在GeoPandas中管理坐标参考系

GeoPandas提供了丰富的函数接口以管理坐标参考系,具体使用如下:

3.3.1 坐标参考系管理

查看坐标管理系

# 框选数据
import geopandas as gpd

# 读取江苏省数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")

# 查看坐标管理系,可以看到默认坐标参考系为EPSG:4326,单位为度。
data.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
# 计算面积,可以看到面积计算结果不对,因为单位是度
# 直接计算面积会报错,提示UserWarning: Geometry is in a geographic CRS
sum(data.area)
9.997823062190006
# 绘图
data.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6e5aed0>

png

设置坐标参考系

对于有些几何形状或者坐标点没有设置坐标参考系,为了告诉GeoPandas如何解释这些几何形状,需要设置默认坐标系。GeoPandas提供set_crs函数来执行这一操作。

p = gpd.GeoSeries([geometry.Point([134.1233,23.456])])
# 没有默认坐标系
p.crs
# 重新设置坐标系
p = p.set_crs("EPSG:4326")
p.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

重设坐标参考系

GeoPandas提供了crs函数以重新设置坐标参考系。

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
world.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6d79d90>

png

# 将坐标参考系转换为EPSG:3857
# 这里舍去北极圈和南极圈数据,在这些地区WGS84/ Pseudo-Mercator变形严重
world = world.to_crs('EPSG:3857')
world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
world.crs
<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World between 85.06°S and 85.06°N.
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
# 横纵坐标尺度变为米
world.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6d65c10>

png

3.3.2 选择合适的投影坐标参考系

坐标参考系选择方法

如果我们获得是地理坐标系数据,如WGS84,在没有更多信息的情况下,如何转换到投影坐标参考系。一个好的办法就是根据所在区域的几何中心点,从前面提到的
espg.io
寻找合适EPSG编码。
这种方法很不专业,只是低精度使用!

以投影江苏省地图为例,一种办法是找到整个国家的几何中心,一种是找到当地的几何中心。

# 框选数据
import geopandas as gpd

# 读取整个国家的数据,无子区域数据
data_china = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/100000.json")

# 查看经纬度几何中心,POINT (103.9 36.4)
# 东经103.9
data_china.centroid
0    POINT (103.9 36.4)
dtype: geometry
# 框选数据
import geopandas as gpd

# 读取江苏省数据,无子区域数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000.json")

# 查看经纬度几何中心,POINT (119.5 33.0)
# 东经119.5
data.centroid
0    POINT (119.5 33.0)
dtype: geometry

接下来我们可以在
espg.io
中,搜索New Beijing kind:PROJCRS。这个搜索词表示搜索以New Beijing为关键词,投影坐标(kind:PROJCRS)的结果。当然以beijing kind:PROJCRS和xian kind:PROJCRS也是可以的。这些都是我国常用的坐标系,不添加kind:PROJCRS也没关系,搜索结果Coordinate reference systems中有相应的选项。然后我们查看地区的几何中心在是否坐标系简介的Area of use中。如下图所示:

比如以整个中国的几何中心(103.9 36.4)查找的结果为
EPSG:4573
,如下所示:

以江苏省的几何中心(119.5 33.0)查找的结果为
EPSG:4586
,如下所示:

坐标参考系转换实例

江苏省的面积为10.72万平方公里,通过EPSG:4573计算的面积为10.84万平方公里,EPSG:4586计算的面积为10.38万平方公里。这两种都是近似估算,实际不会这么粗糙。但是可以看到选择合适的EPSG编码结果都大差不差,想要更准确的结果,可以自行查找其他EPSG编码。

此外,为了展示结果可靠性 我们可以获得某一坐标位置进行绘制。通过地图查询苏州市上的某一经纬度坐标(120.601868 31.332525),如下所示:

# 框选数据
import geopandas as gpd
from shapely import geometry

# 读取江苏省数据
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")


ax = data.plot()
# 绘制点,该坐标点需要预先设置crs
p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.plot(ax=ax, color='red', markersize=100)
<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6c62550>

png

# EPSG:4573下计算面积,EPSG:4573单位尺度为米
# 1平方公里等于10的6次方米,进行计算后结果为10.84万平方公里
data.to_crs("EPSG:4573").area.sum()/1e6/1e4
10.849724500925806
ax = data.to_crs("EPSG:4573").plot()

p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.to_crs("EPSG:4573").plot(ax=ax, color='red', markersize=100)
p.to_crs("EPSG:4573")
0    POINT (19993305.0 3575282.1)
dtype: geometry

png

# EPSG:4586下计算面积
data.to_crs("EPSG:4586").area.sum()/1e6/1e4
10.381446380451157
# 注意绘图结果是有差别的
ax = data.to_crs("EPSG:4586").plot()

p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.to_crs("EPSG:4586").plot(ax=ax, color='red', markersize=100)
p.to_crs("EPSG:4586")
0    POINT (842904.5 3473513.1)
dtype: geometry

png

从上面可以看到EPSG:4573和EPSG:4586表示地图的横坐标单位不一样。EPSG:4573代表的是New Beijing / Gauss-Kruger zone 18,这个zone就是指带号的意思。EPSG:4586代表的是New Beijing / Gauss-Kruger CM 117E,表示没有加带号。简单解释:例如常见的6度分带投影,表示从零度子午线开始,自西向东每个经差6度为一投影带,全球共分60个带,用1,2,3,4,5,... ,60表示。比如没带号的横坐标为500000米。有带号zone 18,表示在没带号的横坐标前加上带号18,变为18500000。具体可以查阅
投影分带

在本例中,如果以江苏省的几何中心查找加带号的坐标系,则应该是New Beijing / Gauss-Kruger zone 20(EPSG:4575)。与EPSG:4586相比,以EPSG:4575为基准绘图,横坐标加了带号(20),面积计算结果和是相等的。这些只是GIS基础知识,本文说得非常简单也有诸多错误(因为作者非专业人士),低精度软件绘制地图也不需要用到过多专业的知识。如果想高精度绘图,请从头开始学习GIS相关知识。

data.to_crs("EPSG:4575").area.sum()/1e6/1e4
10.381446380451155
ax = data.to_crs("EPSG:4575").plot()
p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.to_crs("EPSG:4575").plot(ax=ax, color='red', markersize=100)

# EPSG:4586 为POINT (842904.5 3473513.1)
p.to_crs("EPSG:4575")
0    POINT (20842904.5 3473513.1)
dtype: geometry

png

4 参考

4.1 GeoPandas使用

4.2 相关资料

学习操作系统原理最好的方法是自己写一个简单的操作系统。


在上一讲中我们介绍了屏幕显示的原理,本讲我们来实战一下。

一、向屏幕输出一个字符mbr4.asm

mbr4.asm中的代码如下:

;将屏幕第一行的第一个字符显示为‘G’。
mov ah,0x07 ;黑底白字
mov al,'G'  ;G的ASCII码是0x47,此时ax的值为0x0747。
mov bx,0xb800
mov es,bx
mov [es:0],ax ;文本模式显存地址从0xb8000开始。

stop: ;标号
hlt
jmp stop 

times 510-($-$$) db 0 ;将从上条指令结束到最后2个字节前的空余字节全部置为0。
db 0x55,0xaa

上面代码中的注释比较详细,结合之前介绍过的内容,大家应该能看懂。下面我们来演示一下。
首先我们回顾一下上节课QEMU中默认显示的内容:

从上面截图中可以看到,QEMU中默认显示的第一行第一个字符是‘S’。
下面我们编译运行mbr4.asm。

从上面这个截图可以看到QEMU第一行第一个字符已经变为了字符‘G’,这是我们第一次向屏幕输出字符。

二、将字符显示到屏幕的任意位置mbr5.asm

在默认的文本模式中,一屏能显示25行80列,共2000个字符。每个字符占用2个显存地址,2000个字符共占用4000个显存地址。所以第一屏的显存地址范围是(0xb8000+0)~(0xb8000+4000)。每行显示80个字符,也就是每行对应160个显存地址。在实际使用中,我们需要能将字符输出到屏幕的任意位置。比如上面的例子中,我们将字符输出到了QEMU本来就有字符的地方,这样混在一起不好。我们看到第二行是空白的,我们下面将字符输出到第二行。
mbr5.asm的代码如下:

mov ax,0xb800
mov es,ax
mov ah,0x07
mov al,'G'
mov [es:160],ax ;将字符'G'显示在屏幕第二行第一个字符的位置

stop:
hlt
jmp stop 

times 510-($-$$) db 0
db 0x55,0xaa

编译运行截图如下:

从上面截图可以看到,我们将字符‘G’显示在了QEMU第二行第一个字符的位置。

三、向屏幕显示字符串mbr6.asm

mbr6.asm的代码如下:

mov ax,0xb800
mov es,ax

;在屏幕第2行显示字符串“GrapeOS"
mov ah,0x07 ;ah中的值一直保持不变
mov al,'G'
mov [es:160],ax
mov al,'r'
mov [es:162],ax ;每个字符对应显存中的2个字节,依次递增2个字节。
mov al,'a'
mov [es:164],ax
mov al,'p'
mov [es:166],ax
mov al,'e'
mov [es:168],ax
mov al,'O'
mov [es:170],ax
mov al,'S'
mov [es:172],ax

stop:
hlt
jmp stop 

times 510-($-$$) db 0
db 0x55,0xaa

编译运行截图如下:

从上面截图中可以看到,我们成功的在QEMU屏幕第二行显示出了字符串“GrapeOS”。


本讲视频版地址:
https://www.bilibili.com/video/BV1VY411v7y2/
本教程代码和资料:
https://gitee.com/jackchengyujia/grapeos-course
GrapeOS操作系统QQ群:643474045

之前没接触过Android maui 开发,这几天刚摸索,有些不合理的地方欢迎指出。

首先添加微信sdk的绑定库

nuget 包:Chi.MauiBinding.Android.WeChat

项目地址:https://github.com/realZhangChi/MauiBinding

在项目目录 Platforms---Android 下建立文件夹wxapi,在wxapi文件夹下添加WXEntryActivity.cs,用来接收微信sdk的回调信息(
Android接入指南 | 微信开放文档 (qq.com)
),建好后修改代码为:

namespaceMauiApp7.Platforms.Android.wxapi
{
//将下面的中文"包名"替换成你实际的包名 [Activity(Name = "包名.WXEntryActivity", Exported = true, TaskAffinity = "包名", LaunchMode = LaunchMode.SingleTask, Label = "微信回调")]//特别注意WXEntryActivity 是集成的Activity 不是MauiAppCompatActivity public classWXEntryActivity : Activity, IWXAPIEventHandler
{
private const string APP_ID = "xxxx";//更改为你自己在微信开放平台应用的appId privateIWXAPI api;protected override voidOnCreate(Bundle savedInstanceState)
{
base.OnCreate(savedInstanceState);
api
= WXAPIFactory.CreateWXAPI(this, APP_ID, false);try{
Intent intent
= this.Intent;
api.HandleIntent(intent,
this);
}
catch(Exception e)
{
Console.WriteLine(
"==========执行了OnCreate Exception:" +e.Message);
}
}
public new voidDispose()
{
this.Dispose();
}
protected override voidOnNewIntent(Intent intent)
{
base.OnNewIntent(intent);this.Intent =intent;
api.HandleIntent(intent,
this);

}
public voidOnReq(BaseReq p0)
{


Finish();
}
public voidOnResp(BaseResp p0)
{
//这里接收微信发来的消息 if (p0.Type == 1)
{
if (p0.ErrCode_ == 0)
{
SendAuth.Resp authResp
=(SendAuth.Resp)p0;

}
else{
SendAuth.Resp authResp
=(SendAuth.Resp)p0;

}

}
Finish();

}
}
}

注意:这个WXEntryActivity.cs 中的代码调试时进不去断点,而不是不执行(执行调试时没进断点一直认为不执行里面的方法,后面加了日志打印发现能打印出日志,实际是执行了的),需要手动添加日志打印。这个文件我只是随便加了几个回调信息,它可以接收所有微信的回调,具体的你自己修改。

另外这里接收到微信回调信息后如何通知静态html 我没找到好的解决方法使用了个绕弯路的方法,就是在调需要回调信息的方法时(例如微信登录,支付等),先链接服务器端的websocket,链接后在调用方法,接收到微信回调信息时向远程服务器发送个请求,远程服务器接收到请求通过websocket通知html。如果各位有更好的方法麻烦通知下我谢谢。

在根目录创建一个静态类PublicMethods.cs (类名位置都可以自定义,这个静态类主要给html js 调用使用的,js调用服务端方法
从 ASP.NET Core Blazor 中的 JavaScript 函数调用 .NET 方法 | Microsoft Learn

注意其中的 #if ANDROID  IOS 指在不同的平台下执行操作

usingMicrosoft.JSInterop;usingSystem;usingSystem.Collections.Generic;usingSystem.Diagnostics;usingSystem.Linq;usingSystem.Text;usingSystem.Text.Json;usingSystem.Threading.Tasks;using staticSystem.Runtime.InteropServices.JavaScript.JSType;#if ANDROID
usingAndroid.Content;usingAndroid.Views;usingAndroid.Runtime;usingAndroid.Net;usingCom.Tencent.MM.Opensdk.Modelmsg;usingCom.Tencent.MM.Opensdk.Openapi;usingSystem.Runtime.InteropServices.JavaScript;usingOrg.Json;#elif IOS
usingUIKit;#endif
namespaceMauiApp7
{
public static classPublicMethods
{
     //微信登录
[JSInvokable]
public staticTask WxLogin()
{
#if ANDROIDIWXAPI api= WXAPIFactory.CreateWXAPI(Android.App.Application.Context, "wx00962a9afb38e1e9", true);var isReg = api.RegisterApp("wx00962a9afb38e1e9");if(isReg)
{
SendAuth.Req req
= newSendAuth.Req();
req.Scope
= "snsapi_userinfo"; //只能填 snsapi_userinfo req.State = "wechat_sdk_demo_test";
api.SendReq(req);
}
Console.WriteLine(
"微信登录方法");#endif return Task.FromResult(0);
}
     //微信分享
[JSInvokable]
public staticTask WxShare()
{
#if ANDROIDIWXAPI api= WXAPIFactory.CreateWXAPI(Android.App.Application.Context, "wx00962a9afb38e1e9", true);var isReg = api.RegisterApp("wx00962a9afb38e1e9");if(isReg)
{
string text="test";//初始化一个 WXTextObject 对象,填写分享的文本内容 WXTextObject textObj = newWXTextObject();
textObj.Text
=text;//用 WXTextObject 对象初始化一个 WXMediaMessage 对象 WXMediaMessage msg = newWXMediaMessage();
msg.MediaObject_
=textObj;
msg.Description
=text;

SendMessageToWX.Req req
= newSendMessageToWX.Req();
req.Transaction
= Guid.NewGuid().ToString("N");
req.Message
=msg;
req.Scene
=SendMessageToWX.Req.WXSceneSession;//调用api接口,发送数据到微信 api.SendReq(req);
}
Console.WriteLine(
"微信分享方法");#endif return Task.FromResult(0);
}
}
}

至此服务端代码完毕,下面时html代码

1.新建maui Blazor应用,修改MainPage.xaml 中的代码,删除BlazorWebView 下子内容,修改后的代码为

  <BlazorWebView x:Name="blazorWebView" HostPage="wwwroot/index.html"  BlazorWebViewInitialized="blazorWebView_BlazorWebViewInitialized">
      
    </BlazorWebView>

修改MainPage.xaml.cs中代码添加BlazorWebViewInitialized 事件,此事件是允许BlazorWebView在Android平台下能够同时访问http和https的混合请求,需搭配android:usesCleartextTraffic="true" 使用 具体参考
maui BlazorWebView Android 中混合使用https和http - 落叶子 - 博客园 (cnblogs.com)

  private void blazorWebView_BlazorWebViewInitialized(objectsender, Microsoft.AspNetCore.Components.WebView.BlazorWebViewInitializedEventArgs e)
{
#if ANDROIDe.WebView.Settings.MixedContentMode=Android.Webkit.MixedContentHandling.AlwaysAllow;#endif}

2.将静态html或者打包好的vue、uniapp项目拷贝到wwwroot目录下(先删除此目录下项目生成的文件),因为我测试是用的uniapp打包的html,因此下面都是基于uniapp的方式实现的,其他的都一样的步骤

将uniapp打包后的项目拷贝后修改index.html文件,添加_framework/blazor.webview.js,就是在index.html中的js 文件上面添加blazor.webview.js的引用

 <script src="_framework/blazor.webview.js" autostart="false" crossorigin="anonymous"></script>
    <script src=./static/js/chunk-vendors.22eccfa8.js></script>
    <script src=./static/js/index.5efaff53.js></script>

3.在需要执行微信skd的方法时直接通过js方法调用(微信sdk方法需先在PublicMethods.cs 静态类中注册)
从 ASP.NET Core Blazor 中的 JavaScript 函数调用 .NET 方法 | Microsoft Learn

DotNet.invokeMethodAsync('MauiApp7', 'WxLogin')
.then(data
=>{
console.log(data);
});

这样就可以成功调用微信sdk的方法了

众所周知,ChatGPT可以帮助研发人员编写或者Debug程序代码,但是在执行过程中,ChatGPT会将程序代码的一些相关文字解释和代码段混合着返回,如此,研发人员还需要自己进行编辑和粘贴操作,效率上差强人意,本次我们试图将ChatGPT直接嵌入到代码业务中,让ChatGPT生成可以直接运行的代码。

ChatGPT的主观回答问题

首先,我们向ChatGPT提出一个简单的代码需求:

可以看到,就像上文所描述的那样,ChatGPT会将文字描述和代码片段混合着返回,其实对于真正的需求者来说,文字描述本身是不必要的,因为如果提问者不知道什么是布隆过滤器,也就不会提出布隆过滤器相关的代码需求。

再看ChatGPT返回的布隆过滤器代码:

import hashlib  
  
class BloomFilter:  
    def __init__(self, size, hash_count):  
        self.size = size  
        self.hash_count = hash_count  
        self.bit_array = [False] * size  
  
    def add(self, string):  
        for seed in range(self.hash_count):  
            result = hashlib.sha256((string + str(seed)).encode()).hexdigest()  
            index = int(result, 16) % self.size  
            self.bit_array[index] = True  
  
    def __contains__(self, string):  
        for seed in range(self.hash_count):  
            result = hashlib.sha256((string + str(seed)).encode()).hexdigest()  
            index = int(result, 16) % self.size  
            if not self.bit_array[index]:  
                return False  
        return True

大体上,没有毛病。但是主观性太强,什么是主观性?就是ChatGPT其实不是站在需求者的视角来编写代码,而是站在自己的角度上,它没有考虑业务的上下文关系,也就是类和方法命名、方法参数、以及参数类型或者返回值以及类型,这些东西是否符合需求者当前的代码业务。

当然,这并不是ChatGPT的错,主要是输入的问题描述不够明确和详细,但如果每次都需要将代码业务逻辑转化为文字告诉ChatGPT,又有些画蛇添足,狗尾续貂之感。

基于业务配置ChatGPT

那么怎样将ChatGPT融入业务代码?首先创建Openai接入函数:

import openai  
  
openai.api_key = "apikey"  
  
def generate_code(func, docstring):  
    init_prompt = "You are a Python expert who can implement the given function."  
    definition = f"def {func}"  
    prompt = f"Read this incomplete Python code:\n```python\n{definition}\n```"  
    prompt += "\n"  
    prompt += f"Complete the Python code that follows this instruction: '{docstring}'. Your response must start with code block '```python'."  
  
    response = openai.ChatCompletion.create(  
        model="gpt-3.5-turbo",  
        temperature=0,  
        max_tokens=1024,  
        top_p=1,  
        messages=[  
            {  
                "role": "system",  
                "content": init_prompt,  
            },  
            {  
                "role": "user",  
                "content": prompt,  
            },  
        ],  
    )  
  
    codeblock = response.choices[0].message.content  
    code = next(filter(None, codeblock.split("```python"))).rsplit("```", 1)[0]  
    code = code.strip()  
  
    return code

诀窍就是提前设置好引导词:

init_prompt = "You are a Python expert who can implement the given function."  
    definition = f"def {func}"  
    prompt = f"Read this incomplete Python code:\n```python\n{definition}\n```"  
    prompt += "\n"  
    prompt += f"Complete the Python code that follows this instruction: '{docstring}'. Your response must start with code block '```python'."

这里我们提前设置两个参数func和docstring,也就是函数名和功能描述,要求ChatGPT严格按照参数的输入来返回代码,现在运行函数:

if __name__ == '__main__':  
      
    print(generate_code("test","Sum two numbers"))

程序返回:

➜  chatgpt_write_code /opt/homebrew/bin/python3.10 "/Users/liuyue/wodfan/work/chatgpt_write_code/chatgpt_write_code.p  
y"  
def test(a, b):  
    return a + b

如此一来,ChatGPT就不会返回废话,而是直接交给我们可以运行的代码。

装饰器调用ChatGPT

事实上,函数调用环节也可以省略,我们可以使用Python装饰器的闭包原理,直接将所定义函数的参数和描述传递给ChatGPT,随后再直接运行被装饰的函数,提高效率:

import inspect  
from functools import wraps  
  
def chatgpt_code(func):  
    @wraps(func)  
    def wrapper(*args, **kwargs):  
        signature = f'{func.__name__}({", ".join(inspect.signature(func).parameters)}):'  
        docstring = func.__doc__.strip()  
        code = generate_code(signature, docstring)  
        print(f"generated code:\n```python\n{code}\n```")  
        exec(code)  
        return locals()[func.__name__](*args, **kwargs)  
  
    return wrapper

将方法定义好之后,使用基于ChatGPT的装饰器:

if __name__ == '__main__':  
      
    @chatgpt_code  
    def sum_two(num1,num2):  
        """  
        Sum two numbers.  
        """  
  
    print(sum_two(1,2))

程序返回:

➜  chatgpt_write_code /opt/homebrew/bin/python3.10 "/Users/liuyue/wodfan/work/chatgpt_write_code/chatgpt_write_code.p  
y"  
sum_two(num1, num2):  
generated code:  
  
def sum_two(num1, num2):  
    """  
    Sum two numbers.  
    """  
    return num1 + num2  
  
3

直接将业务逻辑和运行结果全部返回。

那么现在,回到开篇的关于布隆过滤器的问题:

if __name__ == '__main__':  
      
    @chatgpt_code  
    def bloom(target:str,storage:list):  
        """  
        Use a Bloom filter to check if the target is in storage , Just use this func , no more class  
        """  
  
    print(bloom("你好",["你好","Helloworld"]))

程序返回:

➜  chatgpt_write_code /opt/homebrew/bin/python3.10 "/Users/liuyue/wodfan/work/chatgpt_write_code/chatgpt_write_code.p  
y"  
generated code:  
  
def bloom(target, storage):  
    # Initialize the Bloom filter with all zeros  
    bloom_filter = [0] * len(storage)  
      
    # Hash the target and set the corresponding bit in the Bloom filter to 1  
    for i in range(len(storage)):  
        if target in storage[i]:  
            bloom_filter[i] = 1  
      
    # Check if all the bits corresponding to the target are set to 1 in the Bloom filter  
    for i in range(len(storage)):  
        if target in storage[i] and bloom_filter[i] == 0:  
            return False  
      
    return True  
  
True  
➜  chatgpt_write_code

丝滑流畅,和业务衔接得天衣无缝,拉链般重合,不需要挑挑拣拣,也不必复制粘贴。

结语

毫无疑问,ChatGPT确然是神兵利器,吹毛可断,无坚不摧。但工具虽好,也需要看在谁的手里,所谓工具无高下,功力有高深,类比的话,如果倚天剑握在三岁孩童手中,不仅毫无增益,还可能伤其自身,但是握在峨眉掌门灭绝师太手里,那就可以横扫千军如卷席了,那才能体现大宗匠的手段。最后,奉上项目代码,与众乡亲同飨:github.com/zcxey2911/chatgptapi_write_code

k8s容器互联-flannel host-gw原理篇

容器系列文章

容器系列视频

简析host-gw

前面分析了flannel vxlan模式进行容器跨主机通信的原理,但是vxlan模式需要对数据包进行额外的封包解包处理,带来的开销较大。

所以flannel提供了另外一种纯3层转发的通信模式,叫做host-gw,顾明思议,这种模式是将主机作为网关在用了。

先来看下网关在ip通信中的作用,例如,一个tcp包有源ip和目的ip,如果目的ip匹配不到路由信息,那么就会将包转发到网关,在一个发往目的ip的过程中,可能会经过多个网关。

网关的本质是作为ip通信的中转站,网络包在传输过程中,目的ip是不会变的,一直在变化的是mac地址,每到达一台主机,那么目的mac地址就会发生变化,变成下一个网关的mac地址,数据包需要到达的下一台主机被称作”下一跳“(next hop)。

了解了网关的作用,再来看看flannel host-gw模式在k8s节点上做了哪些改动。

集群基本信息

这里我同样是启动了一个3节点的集群,cni插件就是用flannel,模式是host-gw模式。

net-conf.json: |
    {
      "Network": "10.10.0.0/16",
      "Backend": {
        "Type": "host-gw"
      }
    }

集群节点信息

parallels@master:~/k8s$ kubectl get nodes -o wide
NAME      STATUS   ROLES                  AGE   VERSION   INTERNAL-IP    EXTERNAL-IP   OS-IMAGE           KERNEL-VERSION      CONTAINER-RUNTIME
master    Ready    control-plane,master   13d   v1.23.3   192.168.2.17   <none>        Ubuntu 22.04 LTS   5.15.0-58-generic   docker://20.10.12
worker1   Ready    <none>                 13d   v1.23.3   192.168.2.16   <none>        Ubuntu 22.04 LTS   5.15.0-60-generic   docker://20.10.12
worker2   Ready    <none>                 13d   v1.23.3   192.168.2.15   <none>        Ubuntu 22.04 LTS   5.15.0-60-generic   docker://20.10.12

然后用busybox镜像启动了4个pod

parallels@master:~/k8s$ kubectl  get pods -o wide
NAME                       READY   STATUS    RESTARTS   AGE   IP          NODE      NOMINATED NODE   READINESS GATES
busybox-8647b8666c-jpnb6   1/1     Running   0          21m   10.10.1.6   worker1   <none>           <none>
busybox-8647b8666c-pg7ps   1/1     Running   0          21m   10.10.2.4   worker2   <none>           <none>
busybox-8647b8666c-sgf8v   1/1     Running   0          21m   10.10.1.5   worker1   <none>           <none>
busybox-8647b8666c-zlxmm   1/1     Running   0          21m   10.10.2.3   worker2   <none>           <none>

我们的目的就是看看worker1节点上的ip为10.10.1.6 的pod 是如何ping通 worker2节点上的ip为 10.10.2.4 的pod的。

分析集群内部网络流动方向

为了接下来的分析更加形象化,这里我先贴上一张集群内部的网络拓扑图。后续的分析都可以随时回顾下这张图。

网络拓扑图

先从10.10.1.6的pod看起,进入10.10.1.6的pod查看路由信息。

worker1节点上的ip为10.10.1.6的pod路由信息

parallels@master:~/k8s$ kubectl exec -it busybox-8647b8666c-jpnb6 /bin/sh
kubectl exec [POD] [COMMAND] is DEPRECATED and will be removed in a future version. Use kubectl exec [POD] -- [COMMAND] instead.
/ #
/ # ip route
default via 10.10.1.1 dev eth0
10.10.0.0/16 via 10.10.1.1 dev eth0
10.10.1.0/24 dev eth0 scope link  src 10.10.1.6

默认网关是10.10.1.1 ,这个ip地址其实就是worker1节点上cni0网桥的ip地址

可以查到worker1节点上cni0的ip地址

parallels@worker1:~$ ifconfig
cni0: flags=4163<UP,BROADCAST,RUNNING,MULTICAST>  mtu 1500
        inet 10.10.1.1  netmask 255.255.255.0  broadcast 10.10.1.255

所以在ip为10.10.1.6的pod内部去ping上worker2节点的pod ip 10.10.2.4 会匹配上第二条路由信息,然后由eth0网卡出去,网关地址是10.10.1.1,所以网络包就从pod内部传送到了worker1的cni0网桥上。

cni0网桥会将mac地址为其自身mac地址的数据包转发到主机的3层网络中,而具体要怎么路由,则是需要看worker1主机上的路由规则。

parallels@worker1:~$ ip route
default via 192.168.2.1 dev enp0s5 proto dhcp src 192.168.2.16 metric 100
10.10.0.0/24 via 192.168.2.17 dev enp0s5
10.10.1.0/24 dev cni0 proto kernel scope link src 10.10.1.1
10.10.2.0/24 via 192.168.2.15 dev enp0s5
172.17.0.0/16 dev docker0 proto kernel scope link src 172.17.0.1 linkdown
192.168.2.0/24 dev enp0s5 proto kernel scope link src 192.168.2.16 metric 100
192.168.2.1 dev enp0s5 proto dhcp scope link src 192.168.2.16 metric 100

这些节点上路由的配置是由flannel 在每个节点上启动的flanneld进程去进行的配置的,配置信息来源是k8s集群内部的etcd集群

我们发送的数据包目的ip是10.10.2.4 ,它会匹配上worker1主机的第二条路由信息,第二条路由信息是在说访问10.10.0.0/24 网段的数据包都将由enp0s5网卡发出,并且网关地址也就是下一跳的ip地址是192.168.2.17,而192.168.2.17 就是worker2的ip地址。

为了看的更加清晰,我们再来回顾下开局的图。
网络拓扑图

这样数据包就到达到worker2节点了,到了worker2节点后,数据包的如何流动是看worker2节点上的路由规则,所以我们再来看下节点2上面的路由规则。记住数据包的目的ip是10.10.2.4。

parallels@worker2:~$ ip route
default via 192.168.2.1 dev enp0s5 proto dhcp src 192.168.2.15 metric 100
10.10.0.0/24 via 192.168.2.17 dev enp0s5
10.10.1.0/24 via 192.168.2.16 dev enp0s5
10.10.2.0/24 dev cni0 proto kernel scope link src 10.10.2.1
172.17.0.0/16 dev docker0 proto kernel scope link src 172.17.0.1 linkdown
192.168.2.0/24 dev enp0s5 proto kernel scope link src 192.168.2.15 metric 100
192.168.2.1 dev enp0s5 proto dhcp scope link src 192.168.2.15 metric 100

匹配上了第4条路由规则,发往 10.10.2.0/24 的网段的数据包是要被cni0网桥处理的,所以数据包来到了worker2节点上的cni0网桥上,cni0是如何找到要发送的目的ip的veth端口的呢?

pod内部的eth0 网卡其实就是个veth设备,veth设备一端连接在pod的网路命名空间中,一端连接在网桥上,从veth的一端发出去的网络包一定能够被另一端接收。

网桥收到主机发来的数据包后,首先看自身有没有数据包的目的ip的端口记录,如果有,那么就从该端口发送数据包,因为连接的veth设备,所以从端口发送出去后,一定能到达pod的内部,veth设备就像是网线一样。

如果没有记录,那么网桥会向通过arp协议广播帧,得到回应后便能知道端口与ip的映射关系。从而将数据包发往正确的端口。

这样一个数据包就完全的从一台主机通过路由规则到达到了另外一台主机,而主机ip实际上是被当成网关,作为原ip地址的下一跳地址了。

host-gw的优缺点

相比于vxlan模式,因为少了封包解包的操作,会提升数据传输的性能。但由于这是一个纯3层转发的方案,要想主机作为的网关的前提,必须是集群中的两台主机是一个二层连通的环境中。