Basemap简介
Basemap是Python可视化库Matplotlib下的一个工具包,主要功能是绘制二维地图,对于空间数据的可视化非常重要。Basemap本身不会进行任何绘图,但提供了将坐标转换为25个不同地图投影之一的功能。Matplotlib也可以用于绘制变换坐标中的轮廓,图像,向量,线或点。basemap包括GSSH海岸线数据集,以及来自GMT的河流、州和国家边界的数据集。这些数据集可用于在地图上以几种不同的分辨率绘制海岸线,河流和政治边界。basemap底层使用了GeometryEngine-OpenSource(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"; #fix r 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实战:OpenFlights数据可视化
OpenFlights官网提供的都是独立文本,形如airports.dat, airlines.dat, routes.dat,这里使用DataQuest提供了一个打包的.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='offsetpoints', 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='offsetpoints', alpha=0.8, size='x-small') fig
参考链接: