Python地理数据可视化工具Basemap

Basemap简介

Basemap 是 Python 可视化库 Matplotlib 下的一个工具包，主要功能是绘制二维地图，对于空间数据的可视化非常重要。Basemap本身不会进行任何绘图，但提供了将坐标转换为25个不同地图投影之一的功能。 Matplotlib也可以用于绘制变换坐标中的轮廓，图像，向量，线或点。basemap包括GSSH海岸线数据集，以及来自GMT的河流、州和国家边界的数据集。这些数据集可用于在地图上以几种不同的分辨率绘制海岸线，河流和政治边界。basemap底层使用了Geometry Engine-Open Source(GEOS)库，用来将海岸线和边界特征剪切到所需的地图投影区域。此外，basemap还提供读取shapefile的功能。

Basemap的安装

Basemap无法使用pip install basemap进行安装，如果安装了Anaconda则可以通过canda install basemap进行安装。

对于不能安装的情况，可以下载编译好的安装包：https://www.lfd.uci.edu/~gohlke/pythonlibs/。在basemap安装之前需要先安装pyproj。

安装完后执行代码会报如下错误：KeyError: ‘PROJ_LIB’

解决方案：

import os
os.environ["PROJ_LIB"] = r"D:\ProgramData\Anaconda3\pkgs\proj4-5.2.0-h6538335_1006\Library\share"; #fixr
from mpl_toolkits.basemap import Basemap

或者直接在环境变量中进行设置。

Basemap的使用

basemap中的一些实例对象：

contour() : draw contour lines.（画轮廓线）
contourf() : draw filled contours.（画填充后的轮廓线）
imshow() : draw an image.(在地图上画图）
pcolor() : draw a pseudocolor plot.(伪色图）
pcolormesh() : draw a pseudocolor plot (faster version for regular meshes).
plot() : draw lines and/or markers.(在地图上画线绘图）
scatter() : draw points with markers.（在地图上画散点图）
quiver() : draw vectors.（画向量图，三维就是曲面图）
barbs() : draw wind barbs (画风羽图）
drawgreatcircle() : draw a great circle（画大圆航线）

Basemap基本用法

import warnings
warnings.filterwarnings('ignore') 
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap()
map.drawcoastlines()
# plt.show()
plt.savefig('test.png')

代码说明：

第一行代码引入了Basemap库和matplotlib，这两个都是必须的。
这张地图是由Basemap类创建的，这个类包含很多属性。如果使用默认值，地图具有以lat=0和lon=0为中心,使用普通圆柱投影模式（Plate Carrée投影）显示地图。投影公式简单到不能再简单了：x=lon, y=lat。但这个投影既不等角也不等积，特别在高纬度地区，与实际相差很大，所以并不实用。改变地图的投影方式非常简单，只用在Basemap()中加入projection参数和lat_0, lon_0参数。
在这个例子中，用drawcoastlines()方法画出海岸线，海岸线的数据已经默认包含在了库文件中。
最后使用pyplot中的方法显示和保存图片。

更改投影很容易，只下面我们换另外一种投影 ortho，看一下效果。另外我们使用方法 fillcontinents()和 drawmapboundary()来填充海洋和大地的颜色：

map = Basemap(projection='ortho',lat_0=30, lon_0=120)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color= 'coral',lake_color='aqua')
map.drawcoastlines()
map.drawcountries()
plt.show()

Basemap支持的投影类型：

import mpl_toolkits.basemap
print(mpl_toolkits.basemap.supported_projections)

cyl              Cylindrical Equidistant                 
merc             Mercator                                
tmerc            Transverse Mercator                     
omerc            Oblique Mercator                        
mill             Miller Cylindrical                      
gall             Gall Stereographic Cylindrical          
cea              Cylindrical Equal Area                  
lcc              Lambert Conformal                       
laea             Lambert Azimuthal Equal Area            
nplaea           North-Polar Lambert Azimuthal           
splaea           South-Polar Lambert Azimuthal           
eqdc             Equidistant Conic                       
aeqd             Azimuthal Equidistant                   
npaeqd           North-Polar Azimuthal Equidistant       
spaeqd           South-Polar Azimuthal Equidistant       
aea              Albers Equal Area                       
stere            Stereographic                           
npstere          North-Polar Stereographic               
spstere          South-Polar Stereographic               
cass             Cassini-Soldner                         
poly             Polyconic                               
ortho            Orthographic                            
geos             Geostationary                           
nsper            Near-Sided Perspective                  
sinu             Sinusoidal                              
moll             Mollweide                               
hammer           Hammer                                  
robin            Robinson                                
kav7             Kavrayskiy VII                          
eck4             Eckert IV                               
vandg            van der Grinten                         
mbtfpq           McBryde-Thomas Flat-Polar Quartic       
gnom             Gnomonic                                
rotpole          Rotated Pole

具体怎么使用和样式可以自己测试。

稍微复杂一些的地图：

import warnings
warnings.filterwarnings('ignore') 
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

my_map = Basemap(resolution='i', projection='merc', llcrnrlat=10.0, urcrnrlat=55.0, llcrnrlon=60., urcrnrlon=140.0)
my_map.drawcoastlines()
my_map.drawcountries()
my_map.drawmapboundary()
my_map.fillcontinents(color='coral', lake_color='aqua')
my_map.drawmeridians(np.arange(0, 360, 30),labels=[True,True,True,True])
my_map.drawparallels(np.arange(-90, 90, 30),labels=[True,True,True,True])
my_map.drawmapboundary(color='k', linewidth=1.0, fill_color='b', zorder=None, ax=None)
plt.show()

basemap类的resolution参数设置分辨率级别，他们是’c’（原始），’l’（低），’i’（中），’h’（高），’f’（完整）或None（如果没有使用边界）。这个选项很重要，在全局地图上设置高分辨率边界可能非常慢。
我们可以从放大到特定区域，这里的参数是：
- llcrnrlat – 左下角的纬度
- llcrnrlon – 左下角的经度
- urcrnrlat – 右上角的纬度
- urcrnrlon – 右上角的经度
可通过drawmeridians()和drawparallels()绘制经线和纬线。

添加点与线：

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
mymap = Basemap(llcrnrlon=-100.,llcrnrlat=20.,urcrnrlon=20.,urcrnrlat=60.,lat_0=40.,lon_0=-20.,lat_ts=20.)
# nylat, nylon are lat/lon of New York
nylat = 40.78; nylon = -73.98
# lonlat, lonlon are lat/lon of London.
lonlat = 51.53; lonlon = 0.08
mymap.scatter([nylon,lonlon],[nylat,lonlat], linewidths = 10, marker = 'o', color = 'red',s=20, zorder=40)
mymap.drawgreatcircle(nylon,nylat,lonlon,lonlat,linewidth=2,color='b')
mymap.plot([nylon,lonlon],[nylat,lonlat],linewidth=2,color='green',latlon='True')
mymap.drawcoastlines()
mymap.fillcontinents()
mymap.drawparallels(np.arange(10,90,20),labels=[1,1,0,1])
mymap.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
plt.show()

Basemap实战：Open Flights数据可视化

Open Flights官网提供的都是独立文本，形如 airports.dat, airlines.dat, routes.dat，这里使用Data Quest 提供了一个打包的 .db 文件。

import sqlite3
import pandas as pd
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

db = sqlite3.connect('flights.db')

# 获取所有中国机场数据
airports = pd.read_sql_query(
    """
    SELECT DISTINCT cast(latitude AS float) AS lat,
       cast(longitude AS float) lon
    FROM airports
    WHERE country = 'China';
    """,
    db
)

# 获取所有中日航线数据
routes = pd.read_sql_query(
    """
    SELECT DISTINCT cast(sa.latitude AS float) AS src_lat,
       cast(sa.longitude AS float) AS src_lon,
       cast(da.latitude AS float) AS dest_lat,
       cast(da.longitude AS float) AS dest_lon,
       sa.city as src_city,
       da.city as dest_city
    FROM routes
    INNER JOIN airports sa ON sa.id = routes.source_id
    AND sa.country = 'China'
    INNER JOIN airports da ON da.id = routes.dest_id
    AND da.country = 'Japan';
    """,
    db
)

# 加载行政区边界
chn = map.readshapefile('gadm36_CHN_shp/gadm36_CHN_1', 'chn', color='#555555', linewidth=0.1, ax=ax)
jpn = map.readshapefile('gadm36_JPN_shp/gadm36_JPN_1', 'jpn', color='#555555', linewidth=0.1, ax=ax)

# 将国内机场添加到地图
x, y = map(list(airports['lon']),list(airports['lat']))
map.scatter(x, y, 50, marker='.', color='b', alpha=0.1, ax=ax)


# 绘制中日航线
grouped = routes.groupby('src_city')
src_city_list = grouped.groups.keys()

colors_length = np.random.random((len(src_city_list), 3))
colors = dict([(city, color) for city, color in zip(src_city_list, colors_length)])

for key, group in grouped:
    for index, row in group.iterrows():
        map.drawgreatcircle(
            row['src_lon'],
            row['src_lat'],
            row['dest_lon'],
            row['dest_lat'],
            linewidth=2,
            c=colors[key],
            alpha=0.3,
            ax=ax
        )

# 添加出发城市
src_labels = routes.loc[:, ['src_lat', 'src_lon', 'src_city']]
src_labels.drop_duplicates(subset=['src_city'], keep='first', inplace=True)
src_labels.reset_index(drop=True, inplace=True)

Xs, Ys = map(list(src_labels['src_lon']), list(src_labels['src_lat']))

for i, (x, y) in enumerate(zip(Xs, Ys)):
    ax.annotate(str(src_labels['src_city'][i]), (x,y), xytext=(0, 5), textcoords='offset points', alpha=0.8, size='small') 

#添加到达城市
dest_labels = routes.loc[:, ['dest_lat', 'dest_lon', 'dest_city']]
dest_labels.drop_duplicates(subset=['dest_city'], keep='first', inplace=True)
dest_labels.reset_index(drop=True, inplace=True)

Xd, Yd = map(list(dest_labels['dest_lon']), list(dest_labels['dest_lat']))

for i, (x, y) in enumerate(zip(Xd, Yd)):
    ax.annotate(str(dest_labels['dest_city'][i]), (x,y), xytext=(5, 0), textcoords='offset points', alpha=0.8, size='x-small')

fig

参考链接：

Python地理数据可视化工具Basemap

Basemap简介

Basemap的安装

Basemap的使用

Basemap实战：Open Flights数据可视化

Flask学习之Werkzeug

Python模板引擎jinja2

PEP 324：子进程管理 subprocess

发表回复取消回复

Python地理数据可视化工具Basemap

Basemap简介

Basemap的安装

Basemap的使用

Basemap实战：Open Flights数据可视化

相关文章:

Flask学习之Werkzeug

Python模板引擎jinja2

PEP 324：子进程管理 subprocess

发表回复 取消回复

发表回复取消回复