用OSM的AOI裁剪影像的Python自动化方法

OSM（OpenStreetMap）是一个开放的，用户可编辑的世界地图。AOI（Area of Interst）为兴趣面，包含以个区域内的如名称，土地用途等各种信息，在实际应用中我们常常希望得到以AOI划分的要偶感影像。本文将介绍获取OSM的AOI，并进一步获取到AOI区域遥感影像的方法。本文将步骤大致分为获取AOI、获取影像、裁剪影像三步。

代码可以参考项目：https://github.com/psiang/OSMcrop

从OSM获取AOI

获取方式是使用Overpass的API，官网上给了不少示例来参考。在Python中，主要调用方法如下，其中query是请求查询的指令，有特定的格式，可在官网查看参考：

import overpy

# 获取OverPass类
osm_op_api = overpy.Overpass()
# 调用类中的接口
result = osm_op_api.query(query)

即将指令发送给服务器返回查询结果。这个接口调用频率不能过快，且一次获取的Node数目有一定限制。

AOI的获取

从OSM可以获取到三种元素——Node、Way和Relation。我们要得到的AOI属于一种Way，但这个Way需要符合需求。每个元素都有自己的Feature在tags标签下，可以在此基础上删选，例如只保留土地利用为住宅区的区域。更多的Feature可以看这里。

根据官网示例，条件筛选query指令可以在中括号内用“=”后接一个确定的属性，或者“~”后接正则表达式：

[out:csv(number,length)][maxsize:1000000000];
[[:Template:GeocodeArea:Sweden]]->.searchArea;
way["highway"="path"]["foot"!~".*"]["bicycle"!~".*"]["segregated"!~".*"](area.searchArea);
make stat number=count(ways),length=sum(length());
out;

要筛选土地利用中的住宅区可以改成：

way["landuse"="residential"](area.searchArea);

除了筛选Way的属性以外，还要对筛选的地理位置范围进行限制。本文提供以城市名限定范围和矩形限定范围的方式获取两个例子，供大家参考。其中在矩形限定范围的例子中给出了解决接口访问频率和数目限制的方法。

以城市名限定范围

官网参考的城市名限定查询语句如下：

area[name="Česko"];
node[information="guidepost"](area);
out;

Way中的Nodes时用Node的ID表示的，而为了得到Way上每个结点的坐标，需要把Node也获取下来，将范围限制在获取的Way中：

way["landuse"="residential"](area);
node(w);

假设要得到芬兰首都Helsinki中所有住宅区的AOI，结合对Way属性的筛选，我们可以如下构建请求语句的代码。其中key_kist是一个需求属性字典，name是城市名称。当一个关键字下有多个属性要对应时，使用了正则的方式“|”合并起来，否则直接用“=”即可。

# name = Helsinki
# key_list = {"landuse":["residential"]}


query = "[out:json];\n"
query += "area[\"name\"=\"{0}\"]->.searchArea;\n".format(name)
for key in key_list.keys():
    query += "(\n"
    temp = "\"{0}\"".format(key)
    if len(key_list[key]) == 1:
        temp += "=\"{0}\"".format(key_list[key][0])
    else:
        temp += "~\"{0}".format(key_list[key][0])
        for v in range(1, len(key_list[key])):
            temp += "|{0}".format(key_list[key][v])
        temp += "\""
    query += "way[{0}](area.searchArea);\n".format(temp)
    query += "node(w);\n"
    query += ");\n"
query += "out;\n"
osm_op_api = overpy.Overpass()
result = osm_op_api.query(query)

所生成的指令为：

[out:json];
area["name"="Helsinki"]->.searchArea;
(
way["landuse"="residential"](area.searchArea);
node(w);
);
out;

以矩形限定范围

改成矩形只需要将指令中way的圆括号内换成南西北东的纬/经度就可以了。

way["landuse"="residential"](60.1162,24.7522,60.3041,25.2466);

对于矩形限定范围，可以较好地解决接口次数和数目限制问题：

为了解决调用频率太快的问题，一旦查询返回了OverpassTooManyRequests的异常则休眠一段时间重新查询。
为了解决返回结点数目限制地问题，一旦报错结点数目过多，则将当前的矩形范围分成更小的范围来查询，将获取的结果合并即可。这个划分可以分治地进行下去。

def osm_request(lat1, lon1, lat2, lon2):
    flag = True
    while flag:
        try:
            try:
                # 查询
                query = ("[out:json];\n"
                         "(\n"
                         "way[\"landuse\"=\"residential\"]({0},{1},{2},{3});\n"
                         "node(w);\n"
                         ");\n"
                         "out;\n").format(lat1, lon1, lat2, lon2)
                osm_op_api = overpy.Overpass()
                result = osm_op_api.query(query)
            except ConnectionResetError as e:
                # 划分成更小的范围
                print(e)
                result = overpy.Result()
                dim = 3      # 划分成3*3
                for i in range(dim):
                    for j in range(dim):
                        temp = osm_request(lat1 + (lat2 - lat1) / dim * i,
                                            lon1 + (lon2 - lon1) / dim * i,
                                            lat1 + (lat2 - lat1) / dim * (i + 1),
                                            lon1 + (lon2 - lon1) / dim * (i + 1))
                        result.expand(temp)
            flag = False
        except (Exception, overpy.exception.OverpassTooManyRequests) as e:
            # 休眠一段时间重新查询
            print(e)
            time.sleep(60)
    return result

AOI的保存

通过query返回得到的result如果需要获得Node和Way可以直接用以下方式得到：

result.nodes
result.ways

而Way的tag获取方式也类似，注意构成Way的Nodes的获取方式，如果是要id列表必须要越权。获取Node坐标的方式也很简单：

way = result.ways[0]
# 获取Way的feature
tag = way.tags
# 获取Way的ID
id = way.id
# 得到的是像result格式一样的Nodes
nodes = tags.get_nodes()
# 得到的是组成Way的Nodes的ID列表
node_id = tags._node_ids
# 获取Node的纬度和经度
lat, lon = nodes[0].lat, nodes[0].lon

下面以将AOI信息保留至CSV文件为例子。

def way2csv(way, key_list):
    line = ""
    # 只保留需要的信息
    for k in way.tags.keys():
        if k in key_list.keys() and way.tags[k] in key_list[k]:
            line += "%s, %s,%s" % (way.id, k, way.tags[k])
            break
    # 找Node的坐标记录下来，可提前将Node存入字典nodes_list方便查找
    for id in way._node_ids:
        if id in nodes_list:
            line += ",%s,%s" % (nodes_list[id][0], nodes_list[id][1])
    return line


# AOI获取
wayset = result.ways
fway = open("%s_aoi.csv" % name, "w+")
for way in wayset:
    fway.write(way2csv(way, key_list) + "\n")
fway.close()

获取AOI结果展示

AOI获取结果

获取AOI对应遥感影像

参考自项目https://github.com/modyuan/getmap以及博客https://www.cnblogs.com/modyuan/p/12293843.html。

在上一步我们得到了每个AOI上的点的坐标，下面我们希望通过这些坐标获取到AOI对应的遥感影像。根据参考博客，思路可以是通过获取并合并瓦片得到，大致可以分为以下几步：

求最小外包围矩形
将矩形坐标转换为瓦片坐标
获取瓦片
拼接瓦片为影像
导出影像及其信息

求最小外包围矩形

直接遍历AOI的所有点的坐标，找到经纬度最大和最小的即可。

# pos为从csv文件获取到的坐标串，以纬度1、经度1、纬度2、经度2……的顺序排列

# find the bounder of polygon
max_x = pos[1]
min_x = pos[1]
max_y = pos[0]
min_y = pos[0]
for i in range(0, len(pos), 2):
    x = pos[i+1]
    y = pos[i]
    max_x = max(max_x, x)
    max_y = max(max_y, y)
    min_x = min(min_x, x)
    min_y = min(min_y, y)

这样一来（max_x, min_y）是右下角，（min_x, max_y）是左上角。

将矩形坐标转化为瓦片坐标

互联网地图用墨卡托投影生成瓦片。墨卡托投影是地球放在一个圆柱内的投影。瓦片坐标系是将墨卡托投影地图的左上角作为瓦片坐标系原点，右方向为X轴，下方向Y轴。假设瓦片层级为level，瓦片的X轴和Y轴都被分成 $2^{level}$ 段，这样地图上的瓦片总数为 $2^{level}\times2^{level}$ 。更多细节可以参考这里。

墨卡托投影

假设已知WGS84经纬度 $lon$ 和 $lat$ ，需要的瓦片层级为 $z$ ，则经纬度到瓦片的转换公式如下：

$x=\lfloor \frac{lon+180}{360} \times 2^z \rfloor$

$y=\lfloor (\frac{1}{2}-\frac{\ln(\tan(\frac{\pi lat}{180})+\sec(\frac{\pi lat}{180}))}{2\pi})\times 2^z\rfloor$

这个公式可以将左上角和右下角的坐标转换成所需瓦片的左上角和右下角坐标，由此可以得到所需的所有瓦片坐标。

获取瓦片

像谷歌或者百度的瓦片服务器按照一定格式请求即可返回瓦片。谷歌和百度的格式如下，在x、y、z中分别写入瓦片的X轴坐标、Y轴坐标和瓦片层级即可：

"google": "http://mt2.google.cn/vt/lyrs={style}&hl=zh-CN&src=app&x={x}&y={y}&z={z}",
"amap": "http://wprd02.is.autonavi.com/appmaptile?style={style}&x={x}&y={y}&z={z}"

大家可以点击连接查看：
http://mt2.google.cn/vt/lyrs=s&x=214130&y=114212&z=1
http://wprd02.is.autonavi.com/appmaptile?style=6&x=214130&y=114212&z=18

构成这样的URL后，伪装成浏览器发送请求即可，项目原作者在此基础上还做了多线程的下载来加快速度，本文不赘述。这种下载方式只适合小规模的影像，而AOI的区域往往都不大。

HEADERS = {'User-Agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_7_5) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/29.0.1547.76 Safari/537.36'}
header = ur.Request(url,headers=HEADERS)
err=0
while(err<3):
    try:
        data = ur.urlopen(header).read()
    except:
        err+=1
    else:
        return data
raise Exception("Bad network link.")

拼接瓦片为影像

左上角和右下角的瓦片可以得到X轴的瓦片数目 $lenx$ 和Y轴的瓦片数目 $leny$ 。使用Pillow的工具将这些瓦片按照相对位置粘贴到一起。一个瓦片占据256 * 256的大小。

# 新建一个正确大小的画布
outpic = pil.new('RGBA', (lenx * 256, leny * 256))
for i, data in enumerate(datas):
    # 获取到每一个瓦片的数据
    picio = io.BytesIO(data)
    small_pic = pil.open(picio)
    # 将数据写入画布的正确位置
    y, x = i // lenx, i % lenx
    outpic.paste(small_pic, (x * 256, y * 256))

导出影像及其信息

outpic.save(outfile)

导出影像的TIF文件很简单，关键是要导出合并后图像的坐标信息TWF文件方便后续裁剪处理。TWF文件的格式如下，共六行，每行一个数：

表示字母	含义
A	X方向上的象素分辨率
D	X方向的旋转系数
B	Y方向的旋转系数
E	Y方向上的象素分辨率
C	左上角象素中心Y坐标
F	左上角象素中心X坐标

利用TWF文件中的参数可以通过像素坐标 $(x_{pix},y_{pix})$ 推出经纬度坐标，公式如下：

$lon=A\cdot x_{pix}+B\cdot y_{pix}+C \\ lat=D\cdot x_{pix}+E\cdot y_{pix}+F$

TWF公式

为了获取C和F，即左上角的像素坐标，需要将瓦片坐标转回成经纬度坐标，即之前公式的反函数。左上角瓦片坐标转换后的结果即为F和C。

$lon=360\cdot\frac{x}{2^z}-180 \\ lat=\arctan(\sinh(\pi-2\pi\cdot\frac{y}{2^z}))$

为了获取A和E，除了要知道左上角坐标 $(F, C)$ 还需要知道右下角坐标 $(x_{max}, y_{max})$ 。A和E相当于是每个像素增加多少坐标，所以要得到X轴和Y轴的像素总数，用总偏移坐标量除以像素总数即可。

$A=\frac{x_{max}-F}{256\cdot lenx} \\ E=\frac{y_{max}-C}{256\cdot leny}$

至于D和B，由于不涉及到旋转，所以赋值为0即可。

获取影像及其信息展示

按照AOI范围裁剪影像

这一步可以用Pillow自带的接口生成多边形mask，按照这个mask进行裁剪。但是生成多边形的接口需要得到多边形每个点的像素位置，所以先要将获取到的AOI组成的点的坐标转换成像素坐标。

经纬坐标到像素坐标的转换

在刚刚介绍TWF文件格式的时候可以通过像素坐标得到经纬坐标。

$lon=A\cdot x_{pix}+B\cdot y_{pix}+C \\ lat=D\cdot x_{pix}+E\cdot y_{pix}+F$

但是现在需要的是将经纬坐标转换成像素坐标，其实就是解上式的二元方程组：

$x_{pix}=\frac{B(lat-F)-E(lon-C)}{DB-AE} \\ y_{pix}=\frac{A(lat-F)-D(lon-C)}{AE-DB}$

裁剪影像

将AOI的经纬坐标通过上面的转换变成像素坐标后记为polygon。使用ImageDraw可以绘制一个多边形mask，这个mask是二维的，表示每一个像素的状态，在多边形内则像素状态为1，反之为0。只需要将三通道的每一个像素分别乘上mask对应的状态即可完成裁剪。

此处的裁剪方法是裁剪掉的部分像素为(0,0,0)。这里可以针对像素坐标再做一次最小外接矩形，按这个矩形裁剪一次。

# read image as RGB and add alpha (transparency)
im = Image.open(tif_file).convert("RGB")

# convert to numpy (for convenience)
imArray = numpy.asarray(im)

# create mask
maskIm = Image.new('L', (imArray.shape[1], imArray.shape[0]), 0)
ImageDraw.Draw(maskIm).polygon(polygon, outline=1, fill=1)
mask = numpy.array(maskIm)

# assemble new image (uint8: 0-255)
newImArray = numpy.empty(imArray.shape, dtype='uint8')

# colors (three first columns, RGB)
newImArray[:, :, :3] = imArray[:, :, :3]

# crop
newImArray[:, :, 0] *= mask
newImArray[:, :, 1] *= mask
newImArray[:, :, 2] *= mask
newImArray = newImArray[min_y:max_y, min_x:max_x]

# back to Image from numpy
newIm = Image.fromarray(newImArray, "RGB")
newIm.save(to)

裁剪结果展示

裁剪结果