python markever设置 python makepy

       制作专题地图传统上是“适当的”GIS的保留，例如ArcGIS或QGIS。 虽然这些工具使得使用shapefile变得容易，并且暴露了一系列常见的日常GIS操作，但它们并不特别适合探索性数据分析。 简而言之，如果您需要在使用数据制作地图之前获取，重塑和压缩数据，则使用数据分析工具（如Pandas）并将其耦合到绘图库会更容易。 本教程将演示如下使用：

• Pandas
• Matplotlib
• matplotlib Basemap工具包，用于在地图上绘制2D数据
• Fiona，OGR的Python接口
• Shapely（形状），用于分析和操纵平面几何对象
• Descartes（笛卡尔），将所述几何对象转换为matplotlib“补丁”
• PySAL，一个空间分析库

 我在这里使用的方法使用交互式REPL（IPython Notebook）进行数据探索和分析，并使用Descartes包将各个多边形（在本例中为伦敦的病房）渲染为matplotlib补丁，然后将它们添加到matplotlib轴实例中。 我应该强调，许多绘图操作可以更快地完成，但我的目标是演示如何精确控制某些操作，以便实现例如 您想要的精确线宽，颜色，Alpha值或标签位置。

Package installation（包安装）

本教程使用Python 2.7.x，并且需要以下non-stdlib包：

Jupyter
Pandas
Numpy
Matplotlib
Basemap
Shapely
Fiona
Descartes
PySAL

其中一些软件包的安装可能很繁琐，需要大量的第三方依赖项（GDAL＆OGR，C＆FORTRAN77（是的，真的）编译器）。 如果您对Python软件包安装和从源代码构建软件有经验，请在安装所需的软件包之前安装这些依赖项（如果您使用OSX，Homebrew和/或Kyngchaos是有用的，特别是对于GDAL和OGR）。 一个virtualenv，并跳过本节的其余部分。

对于其他人：除了Descartes（笛卡尔）和PySAL之外，Enthought的Canopy（学术用户免费）几乎可以提供您所需的一切。 您可以非常轻松地将它们安装到Canopy用户Python中，有关详细信息，请参阅此支持文章。

运行代码

我发现Jupyter Notebook最适合这样：代码可以在单元格中单独运行，可以轻松纠正错误，图形内联呈现，可以轻松地微调输出。 打开笔记本很简单：从命令行运行jupyter notebook。 这将在您的浏览器中打开一个新笔记本。 Canopy用户：Canopy编辑器将为您完成此操作。

PEP8

我的线路的一些（确切地说是16）是超长的。 我知道。 对不起。

Obtaining a basemap（获取底图）

我们将使用来自Esri Shapefiles的底图，我们将在伦敦地图上绘制数据。 我为此创建了一个shapefile，它在.zip格式下可以在Crown Copyright下获得。 下载它，并将文件解压缩到主项目目录下的名为data的目录中。

Obtaining some data（获取一些数据）

我们将使用相同的数据制作三张地图：伦敦境内的蓝色斑块位置。 为此，我们将从Open Plaques中提供的主XML文件中提取经度，纬度和一些其他功能。 把它拿到这里并把它放在数据目录中。 此文件包含Open Plaques所知道的每个牌匾的数据，但在某些情况下它不完整，并且在我们开始提取有用的子集之前需要清理。

Extracting and cleaning the data（提取和清理数据）

让我们首先导入我们需要的包。 我将根据需要讨论某些库的重要性。

from lxml import etree
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import Normalize
from matplotlib.collections import PatchCollection
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Point, Polygon, MultiPoint, MultiPolygon
from shapely.prepared import prep
from pysal.esda.mapclassify import Natural_Breaks as nb
from descartes import PolygonPatch
import fiona
from itertools import chain

现在，我们将把所有有用的XML数据提取到一个字典中。

tree = etree.parse("data/plaques_20140619.xml")
root = tree.getroot()

output = dict()
output['raw'] = []
output['crs'] = []
output['lon'] = []
output['lat'] = []

for each in root.xpath('/openplaques/plaque/geo'):
    # check what we got back
    output['crs'].append(each.get('reference_system'))
    output['lon'].append(each.get('longitude'))
    output['lat'].append(each.get('latitude'))
    # now go back up to plaque
    r = each.getparent().xpath('inscription/raw')[0]
    if isinstance(r.text, str):
        output['raw'].append(r.text.lstrip().rstrip())
    else:
        output['raw'].append(None)

这将生成一个包含坐标参考系统，经度，纬度和每个斑块记录描述的字典。 接下来，我们将创建一个Pandas DataFrame，删除所有不包含描述的记录，并将long和lat值从字符串转换为浮点数。

df = pd.DataFrame(output)
df = df.replace({'raw': 0}, None)
df = df.dropna()
df[['lon', 'lat']] = df[['lon', 'lat']].astype(float)

现在，我们将打开我们的shapefile，并从中获取一些数据，以便设置我们的底图：

shp = fiona.open('data/london_wards.shp')
bds = shp.bounds
shp.close()
extra = 0.01
ll = (bds[0], bds[1])
ur = (bds[2], bds[3])
coords = list(chain(ll, ur))
w, h = coords[2] - coords[0], coords[3] - coords[1]

我们在这做了两件事：

1.提取了地图边界
2.计算我们的底图的范围，宽度和高度

我们已经准备好创建一个Basemap实例，我们可以使用它来绘制我们的地图。

m = Basemap(
    projection='tmerc',
    lon_0=-2.,
    lat_0=49.,
    ellps = 'WGS84',
    llcrnrlon=coords[0] - extra * w,
    llcrnrlat=coords[1] - extra + 0.01 * h,
    urcrnrlon=coords[2] + extra * w,
    urcrnrlat=coords[3] + extra + 0.01 * h,
    lat_ts=0,
    resolution='i',
    suppress_ticks=True)
m.readshapefile(
    'data/london_wards',
    'london',
    color='none',
    zorder=2)

我选择了横向墨卡托投影，因为它在东西向范围较小的区域表现出较小的失真。 这个投影要求我们指定一个中心经度和纬度，我将其设置为-2,49。

# set up a map dataframe
df_map = pd.DataFrame({
    'poly': [Polygon(xy) for xy in m.london],
    'ward_name': [ward['NAME'] for ward in m.london_info]})
df_map['area_m'] = df_map['poly'].map(lambda x: x.area)
df_map['area_km'] = df_map['area_m'] / 100000

# Create Point objects in map coordinates from dataframe lon and lat values
map_points = pd.Series(
    [Point(m(mapped_x, mapped_y)) for mapped_x, mapped_y in zip(df['lon'], df['lat'])])
plaque_points = MultiPoint(list(map_points.values))
wards_polygon = prep(MultiPolygon(list(df_map['poly'].values)))
# calculate points that fall within the London boundary
ldn_points = filter(wards_polygon.contains, plaque_points)

请注意，map_points系列是通过将经度和纬度值传递给我们的Basemap实例m来创建的。 这会将坐标从长度和纬度转换为地图投影坐标（以米为单位）。 我们的df_map数据框现在包含以下列：

shapefile中每个病房的多边形
它的描述
面积以平方米计
它的面积是平方公里

我们还从组合的病房多边形中创建了一个准备好的几何对象。 我们这样做是为了加快我们的会员检查操作。 我们通过从map_points创建MultiPolygon，然后使用contains（）方法进行过滤来执行成员资格检查，该方法是返回wards_polygon中包含的所有点的二元谓词。 结果是一个Pandas系列，ldn_points，我们将用它来制作我们的地图。

下面的两个函数可以更轻松地为我们的地图生成彩条。 详细了解文档字符串 - 实质上，其中一个对颜色渐变进行了分类，而另一个则更容易标记颜色条。

# Convenience functions for working with colour ramps and bars
def colorbar_index(ncolors, cmap, labels=None, **kwargs):
    """
    This is a convenience function to stop you making off-by-one errors
    Takes a standard colour ramp, and discretizes it,
    then draws a colour bar with correctly aligned labels
    """
    cmap = cmap_discretize(cmap, ncolors)
    mappable = cm.ScalarMappable(cmap=cmap)
    mappable.set_array([])
    mappable.set_clim(-0.5, ncolors+0.5)
    colorbar = plt.colorbar(mappable, **kwargs)
    colorbar.set_ticks(np.linspace(0, ncolors, ncolors))
    colorbar.set_ticklabels(range(ncolors))
    if labels:
        colorbar.set_ticklabels(labels)
    return colorbar

def cmap_discretize(cmap, N):
    """
    Return a discrete colormap from the continuous colormap cmap.

        cmap: colormap instance, eg. cm.jet. 
        N: number of colors.

    Example
        x = resize(arange(100), (5,100))
        djet = cmap_discretize(cm.jet, 5)
        imshow(x, cmap=djet)

    """
    if type(cmap) == str:
        cmap = get_cmap(cmap)
    colors_i = np.concatenate((np.linspace(0, 1., N), (0., 0., 0., 0.)))
    colors_rgba = cmap(colors_i)
    indices = np.linspace(0, 1., N + 1)
    cdict = {}
    for ki, key in enumerate(('red', 'green', 'blue')):
        cdict[key] = [(indices[i], colors_rgba[i - 1, ki], colors_rgba[i, ki]) for i in xrange(N + 1)]
    return matplotlib.colors.LinearSegmentedColormap(cmap.name + "_%d" % N, cdict, 1024)

让我们做一个散点图

# draw ward patches from polygons
df_map['patches'] = df_map['poly'].map(lambda x: PolygonPatch(
    x,
    fc='#555555',
    ec='#787878', lw=.25, alpha=.9,
    zorder=4))

plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, axisbg='w', frame_on=False)

# we don't need to pass points to m() because we calculated using map_points and shapefile polygons
dev = m.scatter(
    [geom.x for geom in ldn_points],
    [geom.y for geom in ldn_points],
    5, marker='o', lw=.25,
    facecolor='#33ccff', edgecolor='w',
    alpha=0.9, antialiased=True,
    label='Blue Plaque Locations', zorder=3)
# plot boroughs by adding the PatchCollection to the axes instance
ax.add_collection(PatchCollection(df_map['patches'].values, match_original=True))
# copyright and source data info
smallprint = ax.text(
    1.03, 0,
    'Total points: %s\nContains Ordnance Survey data\n$\copyright$ Crown copyright and database right 2013\nPlaque data from http://openplaques.org' % len(ldn_points),
    ha='right', va='bottom',
    size=4,
    color='#555555',
    transform=ax.transAxes)

# Draw a map scale
m.drawmapscale(
    coords[0] + 0.08, coords[1] + 0.015,
    coords[0], coords[1],
    10.,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#555555',
    fontcolor='#555555',
    zorder=5)
plt.title("Blue Plaque Locations, London")
plt.tight_layout()
# this will set the image width to 722px at 100dpi
fig.set_size_inches(7.22, 5.25)
plt.savefig('data/london_plaques.png', dpi=100, alpha=True)
plt.show()

我们在地图上绘制了一个散点图，其中包含直径为50米的点，对应于数据框中的每个点。

这是第一步，但并没有真正告诉我们每个病房密度的任何有趣之处 - 仅仅是在伦敦市中心发现的斑块多于外围病房。

Creating a Choropleth Map, Normalised by Ward Area（创建一个由病区标准化的等值区域图）

df_map['count'] = df_map['poly'].map(lambda x: int(len(filter(prep(x).contains, ldn_points))))
df_map['density_m'] = df_map['count'] / df_map['area_m']
df_map['density_km'] = df_map['count'] / df_map['area_km']
# it's easier to work with NaN values when classifying
df_map.replace(to_replace={'density_m': {0: np.nan}, 'density_km': {0: np.nan}}, inplace=True)

我们现在已经创建了一些额外的列，其中包含每个病房中的点数，以及每个病房的每平方米和平方公里的密度。 像这样的标准化允许我们比较病房。

我们几乎准备好制作一个等值线图，但首先，我们必须将我们的病房分成几类，以便轻松区分它们。 我们将使用名为Jenks Natural Breaks的迭代方法来实现这一目标。

# Calculate Jenks natural breaks for density
breaks = nb(
    df_map[df_map['density_km'].notnull()].density_km.values,
    initial=300,
    k=5)
# the notnull method lets us match indices when joining
jb = pd.DataFrame({'jenks_bins': breaks.yb}, index=df_map[df_map['density_km'].notnull()].index)
df_map = df_map.join(jb)
df_map.jenks_bins.fillna(-1, inplace=True)

我们计算了包含一个或多个斑块（density_km不是Null）的所有病房的类（在本例中为五个），并创建了一个包含类号（0  -  4）的新数据帧，其索引与 非零密度值。 这样可以轻松地将其加入现有数据框。 最后一步涉及将bin类-1分配给所有非值行（wards），以便创建单独的零密度类。

我们还想为我们的课程创建一个合理的标签：

jenks_labels = ["<= %0.1f/km$^2$(%s wards)" % (b, c) for b, c in zip(
    breaks.bins, breaks.counts)]
jenks_labels.insert(0, 'No plaques (%s wards)' % len(df_map[df_map['density_km'].isnull()]))

这将显示密度/平方公里，以及班级中的病房数量。

我们现在准备绘制我们的等值区域图：

plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, axisbg='w', frame_on=False)

# use a blue colour ramp - we'll be converting it to a map using cmap()
cmap = plt.get_cmap('Blues')
# draw wards with grey outlines
df_map['patches'] = df_map['poly'].map(lambda x: PolygonPatch(x, ec='#555555', lw=.2, alpha=1., zorder=4))
pc = PatchCollection(df_map['patches'], match_original=True)
# impose our colour map onto the patch collection
norm = Normalize()
pc.set_facecolor(cmap(norm(df_map['jenks_bins'].values)))
ax.add_collection(pc)

# Add a colour bar
cb = colorbar_index(ncolors=len(jenks_labels), cmap=cmap, shrink=0.5, labels=jenks_labels)
cb.ax.tick_params(labelsize=6)

# Show highest densities, in descending order
highest = '\n'.join(
    value[1] for _, value in df_map[(df_map['jenks_bins'] == 4)][:10].sort().iterrows())
highest = 'Most Dense Wards:\n\n' + highest
# Subtraction is necessary for precise y coordinate alignment
details = cb.ax.text(
    -1., 0 - 0.007,
    highest,
    ha='right', va='bottom',
    size=5,
    color='#555555')

# Bin method, copyright and source data info
smallprint = ax.text(
    1.03, 0,
    'Classification method: natural breaks\nContains Ordnance Survey data\n$\copyright$ Crown copyright and database right 2013\nPlaque data from http://openplaques.org',
    ha='right', va='bottom',
    size=4,
    color='#555555',
    transform=ax.transAxes)

# Draw a map scale
m.drawmapscale(
    coords[0] + 0.08, coords[1] + 0.015,
    coords[0], coords[1],
    10.,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#555555',
    fontcolor='#555555',
    zorder=5)
# this will set the image width to 722px at 100dpi
plt.tight_layout()
fig.set_size_inches(7.22, 5.25)
plt.savefig('data/london_plaques.png', dpi=100, alpha=True)
plt.show()

最后，我们可以使用hex bin创建一个替代地图。 正如我们将要看到的，这些是点图的更具信息性的替代方案。 Basemap包提供了一个hex-binning方法，我们需要一些额外的信息才能使用它：

1.将要使用的点的经度和纬度坐标必须作为numpy数组提供。

2.我们必须指定网格大小，以米为单位。 您可以尝试此设置; 我选择了125。

3.将mincnt值设置为1意味着在网格中没有斑块的区域中不会绘制任何分档。

4.您可以指定bin类型。 在这种情况下，我选择了log，它使用颜色映射的对数标度。 这更清楚地强调了每个箱子的密度的微小差异。

The code:

# draw ward patches from polygons
df_map['patches'] = df_map['poly'].map(lambda x: PolygonPatch(
    x, fc='#555555', ec='#787878', lw=.25, alpha=.9, zorder=0))

plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, axisbg='w', frame_on=False)

# plot boroughs by adding the PatchCollection to the axes instance
ax.add_collection(PatchCollection(df_map['patches'].values, match_original=True))

df_london = df[
    (df['lon'] >= ll[0]) &
    (df['lon'] <= ur[0]) &
    (df['lat'] >= ll[1]) &
    (df['lat'] <= ur[1])]

lon_ldn = df_london.lon.values
lat_ldn = df_london.lat.values

# the mincnt argument only shows cells with a value >= 1
# hexbin wants np arrays, not plain lists
hx = m.hexbin(
    np.array([geom.x for geom in ldn_points]),
    np.array([geom.y for geom in ldn_points]),
    gridsize=125,
    bins='log',
    mincnt=1,
    edgecolor='none',
    alpha=1.,
    lw=0.2,
    cmap=plt.get_cmap('Blues'))

# copyright and source data info
smallprint = ax.text(
    1.03, 0,
    'Total points: %s\nContains Ordnance Survey data\n$\copyright$ Crown copyright and database right 2013\nPlaque data from http://openplaques.org' % len(ldn_points),
    ha='right', va='bottom',
    size=4,
    color='#555555',
    transform=ax.transAxes)

# Draw a map scale
m.drawmapscale(
    coords[0] + 0.08, coords[1] + 0.015,
    coords[0], coords[1],
    10.,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#555555',
    fontcolor='#555555',
    zorder=5)

plt.title("Blue Plaque Density, London")
plt.tight_layout()
# this will set the image width to 722px at 100dpi
fig.set_size_inches(7.22, 5.25)
plt.savefig('data/london_plaques.png', dpi=100, alpha=True)
plt.show()

您可以在此处使用IPython Notebook Viewer查看和下载代码的工作副本。 请注意，您必须有一个名为data的子目录，包含XML数据源和Shapefile才能运行它。

在以后的文章中，我将讨论Geopandas ......