python安装basemap

basemap工具集（https://matplotlib.org/basemap/）使得我们能够用python在地图上绘制2D数据。basemap提供了许多不同的地球投影一级一种将地球上的经纬度坐标投影转换为二维matplotlib图的方式

1. 打开网址（点这里），下载basemap和pyproj这两个安装包，找到适应自己电脑和python的对应版本
2. 安装

​ 先执行pip install pyproj-2.6.0-cp36-cp36m-win_amd64.whl

​ 再执行 pip install basemap-1.2.1-cp36-cp36m-win_amd64.whl

3. 如果看到提示成功，说明安装完成

绘制地图仪截面

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

#绘制地理信息相关的图像，当参数 projection 的值为 'ortho'时，将得到一个地球仪截面
plt.figure(figsize=(8, 8))
m = Basemap(projection='ortho', resolution=None, lat_0=50, lon_0=-100)
m.bluemarble(scale=0.5)
plt.show()

运行结果

图形化显示海地地震危机数据

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

# 绘制地图：图形化显示海地地震危机数据
path =str(os.getcwd()).replace('\\','/')
filename = path+'/datasets/ch08_data/Haiti.csv'
data = pd.read_csv(filename,engine='python') # 2010年海地地震及其余震期间搜集的数据
# 清除错误位置信息并移除缺失分类信息
data = data[(data.LATITUDE > 18) & (data.LATITUDE <20) & (data.LONGITUDE>-75) & (data.LONGITUDE<-70 ) & data.CATEGORY.notnull()]

# 获取所有分类的列表
def get_all_categories(cat_series):
    cat_sets = [set(to_cat_list(x)) for x in cat_series]
    return sorted(set.union(*cat_sets))

def to_cat_list(catstr):
    stripped = (x.strip() for x in catstr.split(',')) # 某一条数据按逗号拆分开来
    return [x for x in stripped if x]

# 将各个分类信息拆分为编码和英语名称
def get_english(cat):
    code,names = cat.split('.')
    if '|' in names:
        names = names.split('|')[1]
    return code,names.strip()

# 将编码和名称映射起来的字典
all_cats = get_all_categories(data.CATEGORY)
english_mapping = dict(get_english(x) for x in all_cats)

# 构造一个全0的DataFrame
def get_code(seq):
    return [x.split('.')[0] for x in seq if x ]
all_codes = get_code(all_cats)
code_index = pd.Index(np.unique(all_codes))
dummy_frame = pd.DataFrame(np.zeros((len(data),len(code_index))),index = data.index,columns=code_index)

# 将各行中适当的项设置为1，再与data进行连接
for row, cat in zip(data.index,data.CATEGORY):
    codes = get_code(to_cat_list(cat))
    dummy_frame.ix[row,codes]=1
data = data.join(dummy_frame.add_prefix('category_'))

def basic_haiti_map(ax=None,lllat=17.25,urlat=20.25,lllon=-75,urlon=71):
    # 创建极球面投影的Basemap实例
    m = Basemap(ax=ax,projection='stere',
                lon_0=(urlon+lllon)/2,
                lat_0=(urlat+lllat)/2,
                llcrnrlat=lllat,urcrnrlat=urlat,
                llcrnrlon=lllon,urcrnrlon=urlon,
                resolution='f')
    
    # 绘制海岸线，州界，国界以及地图边界
    m.drawcoastlines()
    m.drawstates()
    m.drawcountries()
    return m

# 让返回的basemap对象知道将坐标转换到画布上
fig,axes = plt.subplots(nrows=2,ncols=2,figsize=(12,10))
fig.subplots_adjust(hspace=0.05,wspace=0.05)

to_plot=['2a','1','3c','7a']
lllat=17.25
urlat=20.25
lllon=-75
urlon=-71

for code,ax in zip(to_plot,axes.flat):
    m = basic_haiti_map(ax,lllat=lllat,urlat=urlat,lllon=lllon,urlon=urlon)
    cat_data = data[data['category_%s' % code]==1]
    
    # 计算地图的投影坐标
    x,y = m(cat_data.LONGITUDE.values,cat_data.LATITUDE.values)
    
    m.plot(x,y,'k.',alpha=0.5)
    ax.set_title('%s: %s' % (code,english_mapping[code]))运行

运行结果

从上图中发现，大部分数据都集中在人口稠密的城市–太子港。basemap还可以叠加来自shapefile的地图数据。首先需要下载一个带有太子港道路的shapefile，然后运行下面脚本即可得到反映事物短缺情况的图片

fig,ax = plt.subplots(nrows=1,ncols=1,figsize=(12,10))
fig.subplots_adjust(hspace=0.05,wspace=0.05)

lllat = 18.43;urlat = 18.69;lllon = -72.57;urlon = -72.08
m = basic_haiti_map(ax,lllat=lllat,urlat=urlat,lllon=lllon,urlon=urlon)

shapefile_path=path+'/datasets/ch08_data/PortAuPrince_Roads/PortAuPrince_Roads'
m.readshapefile(shapefile_path,'road')

code='2a'
cat_data = data[data['category_%s' % code]==1]
# 计算地图的投影坐标
x,y = m(cat_data.LONGITUDE.values,cat_data.LATITUDE.values)
m.plot(x,y,'k.',alpha=0.5)

运行结果