roof detection by opencv

coding
amap api, aws, google map api, opencv, python

背景简介

可再生能源是世界的未来，光伏发电是可再生能源里不可或缺的一个版块。
本人有幸自14年从传统石化能源行业转行到新能源行业，并且一直在与分布式光伏项目打交道。

对于分布式项目国内外也有不同，北美的分布式项目多指容量较小的项目，也包含部分直接对电网ppa的项目。
而我国的分布式项目主要是指屋顶分布式项目，虽然走量大的是全额上网项目，但是个人还是坚持自发自用项目才是真正的分布式项目，
也是最符合业务逻辑(高风险高回报)的项目。

那么好的屋顶(面积要大，承重要够)和好的企业(讲信用，按时付电费，存续时间长)是我们主要搜索的目标。

因此有足够的动机通过卫星图挖掘好的屋顶资源，抢占先机。虽说拿下屋顶主要靠人脉资源，
但是发掘优质资源勉强算是一种需求。

该应用主要用于实现下列事项：

  1. 下载选定区域卫星图
  2. 识别图中的屋顶
  3. 获取屋顶的业主企业信息

可用的免费卫星图源也不多，就google, 百度，高德，其他付费图源包括HERE, Openstreet map以及其他国内付费卫星图源。
就时效性和可用性(API配额、便捷和技术支持)而言就只剩google了。唯一麻烦的是需要翻墙，在自己电脑上使用shadowsocks
全局下载很不稳定，1600张600x650 z16的图断断续续地下了将近一周。

识别程序主要参考了 Roof detection in aerial images of Uganda, 该文献google暂时查不到，scholar里也没有，比较奇怪。
该文献内容虽然比较粗糙，但是讲的方法是最直白最实用的了。

Roof detection in aerial images of Uganda,
Angela Santin Ceballos, Master of Science,
Artificial Intelligence School of Informatics,
University of Edinburgh, 2015

文献中功能实现的代码详见github. 其pipeline详见下文：

1
2
3
4
5
6
opencv cascading training -> cascade models : fast propose rect roof candidates, output the coords of the roofs in satellite image
python machine learning pkg Theano -> convolutional neural network model : recognite whether it's a roof including rotated scenario in candidates, output the possibility
raw image ->> opencv viola jones detection with cascede models -> cNN detection -> Non-maxima suppression ->> detected image  
Non-maxima suppression used to delete the roof parts and save the most overfeat one

获取企业信息主要使用map api服务商的reverse geocoding功能，即逆地理信息查询，通过经纬度查询point of interest(poi).
poi的结果严重依赖地图服务商的数据库。这方面就是国内地图服务商的强项了，google显然离开中国市场太久导致无人上报，水土不服。
后经过斟酌选用高德API来进行reGeo,主要原因是配额和技术支持以及地点的可靠性。高德个人开发者每天只有2000配额，
升级为企业开发者后拥有每天400万配额，每分钟6万上限。缺点是识别出来的图片的坐标系统为WSG89,且国内有偏移，
需要用高德的坐标准换API转换后才能减小偏移误差，提升poi匹配精度。即使这样也无法完全消除坐标偏移影响。

由于该项目初期对ap要求不高，因此只实现了前面opencv cascading detection部分，后续cNN未实施。
初步计划在本地执行python脚本,
但是在实际过程中发现脚本的运算时间过长(dell xps13 9350, ubuntu 16 LTS 平均识别时间约20s),效率太低。
后换闲置的dell vostro 3800挂机运行，仅执行识别脚本，每天最多只能输出~9000 poi, ~1200 km2硬是跑了5天多，效率依旧很低。

后使用aws ec2 t2.micro配置的服务器进行图像下载和识别，大概每小时可以处理约100 km2区域，涵盖约3000 poi.

代码分析

卫星图像获取部分代码详见下文：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
lat_list = np.arange(start_lat, end_lat, step)
lon_list = np.arange(start_lon, end_lon, step)
total_num = len(lat_list) * len(lon_list)
a = np.array([[0, 0]], dtype=float)
for lat in np.arange(start_lat, end_lat, step):
    for lon in np.arange(start_lon, end_lon, step):
        b = np.array([[lat, lon]], dtype=float)
        a = np.concatenate((a, b))
flag = 2
i = flag
for v in a[flag:]:
    lat = v[0]
    lon = v[1]
    urlparams = urllib.urlencode({'center': '{0},{1}'.format(lat, lon),
                                  'zoom': 16,
                                  'size': '600x625',
                                  'maptype': 'satellite',
                                  'key': GOGL_WEB_KEY})
    url = 'https://maps.googleapis.com/maps/api/staticmap?' + urlparams
    print 'now saving... {0}/{1} ...'.format(i, total_num)
    urllib.urlretrieve(url, "image_cache/sate_c_{0}_{1}_z_16.png".format(lat, lon))
    with open('scan_coords.csv', 'ab') as f:
        writer = csv.writer(f, delimiter=';')
        writer.writerow([i, "image_cache/sate_c_{0}_{1}_z_16.png".format(lat, lon), lat, lon])
    i += 1
    t = random.randint(1, 3)
    time.sleep(t)
print '*** mission complete ***'

其中的难点主要在于经纬度入参，需要根据方块边界的经纬度每张图片的经纬度中心点，
下载像素范围(与zoom有关)以及步长。经纬度包含小数点后6位，需要进行float除法，并且组成数组遍历。
但是list只支持int, 因此需要使用numpy中的array来完成按步长计算和遍历。

1
2
3
4
5
6
7
8
input_xml = 'parameter/cascade.xml'
roof_cascade = cv2.CascadeClassifier(input_xml)
img = cv2.imread(input_image)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
roofs = roof_cascade.detectMultiScale(gray, 1.05, 3)

其中Input xml即为训练完成的model,识别代码也非常简单，roofs为输出的obj, 内含像素坐标集合。

cascading训练相关资料

在初次训练的时候花费了近2天的时间手工截取了1200多个正样本，训练效果一般。这种做法比较不可取，但是也局限于没法找到免费的训练样本集。
上文中的作者直接pitch了在Uganda大学做相关研究的教授，直接获取了训练样本，节省了大量的时间。

object detection case及相关资料

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def regeo(lat, lon):
    urlparams = urllib.urlencode({'key': AMAP_WEB_KEY,
                                  'location': '{0},{1}'.format(lon, lat),
                                  'poitype': '公司',
                                  'radius': 300,
                                  'extensions': 'all',
                                  'batch': 'false',
                                  'roadlevel': 1})
    url = 'http://restapi.amap.com/v3/geocode/regeo?' + urlparams
    urlhandler = urllib.urlopen(url)
    mdata = json.loads(urlhandler.read())
    addr = mdata['regeocode']['formatted_address']
    company = list()
    company.append(addr)
    note = 'there does not exist any companies near 300m of this coords'
    for poi in mdata['regeocode']['pois']:
        if u'公司' in poi['type']:
            company.append(poi['name'])
    if company:
        return company
    else:
        return note
        print note

这里只是调用高德的reGeo API, 并没有什么特别要注意的难点，内里也只有一些简单的string manipulate.

trouble shooting

图像识别中处理旋转的图像和识别准确率的提升

后续程序的优化在于使用cNN方法进行图像识别，opencv的VJ方法的有点是速度快，可以快速推送candidates,
后跟进cNN进行屋顶的是非概率判断也是非常关键的。当然在VJ方法中也可以选装图形进行训练或者识别，
毕竟屋顶的朝向多数为正南，东南或者西南，特征也非常显著。总之存在很多方法提升ap.

也是自己挖的坑还需要再填补。瓶颈在于训练样本的获取。

识别效率

现阶段看来使用aws ec2的服务可以完全满足需求，商用价格也相对合适，自用的话第一次使用有1年的免费t2.micro的使用期。
付费的可以选择light sail服务，$10月租与t2.micro同配置的服务器也算合适。

国内的BAT三家也都提供相应的云服务，由于使用阿里旗下的高德地图API, 因此尝试过使用阿里云的计算资源，后由于没有免费的试用资源放弃。
aws中国是单独代理运营的公司，且不提供试用，要备案，约束较多，不建议使用。腾讯的云服务有争议，百度的就不推荐使用了。

分享到 评论