基于百度地图API在AI Studio上的卫星地图块图像处理与分类

本项目基于百度地图API获取了不同的瓦片地图并进行合并等处理，可用于遥感和抽象地图的地图块的图像分类、分割、检测等数据的制作。

AI Studio

920人浏览 · 2022-12-28 22:53:27

AI Studio · 2022-12-28 22:53:27 发布

基于百度地图API在AI Stduio上的瓦片地图块图像处理与分类

项目介绍

本项目基于百度地图API获取了不同的瓦片地图并进行合并等处理，可用于遥感和抽象地图的地图块的图像分类、分割、检测等数据的制作。

项目背景

卫星地图的图像分类是依据是地物的光谱特征，确定判别函数和相应的判别准则，将图像所有的像元按性质分为若干类别的过程，主要方式分为监督分类与非监督分类。
深度学习( deep learning) 作为机器学习算法中的一个新兴技术，其动机在于建立模拟人脑进行分析学习的神经网络，它能通过海量的训练数据和具有很多隐藏层的深度模型学习更有用的特征，最终提升分类的准确性。近年来深度学习在图像分类应用中取得了令人瞩目的成绩，越来越多的学者开始将深度学习应用于遥感图像处理中。几种常用的深度学习方法包括自动编码器、卷积神经网络、深度信念网络和针对小训练样本的迁移学习。这里讨论的方法即为使用百度地图API获取卫星地图，从而进行相应的图像分类与应用。

百度地图开放平台

申请API

1、打开网页 http://lbsyun.baidu.com/index.php?title=首页

2、使用之前需要注册成为开发者，根据提示进行实名认证，注册成功即可。

3、注册成功之后，打开控制台，进入应用管理，输入应用名称，类型，白名单写英文半角的*。

4、创建成功，显示下图，会显示密钥AK

简单测试验证-1

根据上面的操作后，设置AK可进行简单的验证。

(说明：需要修改经度和维度可从百度地图的坐标拾取系统进行获取，
链接为：https://api.map.baidu.com/lbsapi/getpoint/index.html)

import requests
import json

#'复制上文中的AK'
ak = 'dYg9wxs2v2xfaNEpTmcR4S6jZ0idkHN6'

origin = '40.01116,116.339303'  # 起始点经纬度坐标
destination = '39.936404,116.452562'  # 终点经纬度坐标
url = 'https://api.map.baidu.com/directionlite/v1/walking?origin=' \
      + origin + '&destination=' + destination + '&ak=' + ak + '&coord_type=wgs84'

'''
url数据解释
1. walking是步行导航的意思，如果想要其他导航方式可参考官方文档
   riding是骑行导航
   driving是驾车导航
   transit是公交导航等等
2. origin是起始点的经纬度坐标，前者为纬度，后者为经度
3. destination是终点的经纬度坐标
4. ak是就是上文中的ak
5. coord_type=wgs84的意思是gps坐标。因为默认是百度经纬度坐标。
'''

response = requests.get(url)  # 请求连接
answer = response.json()  # 获取json格式
answer_json = json.dumps(answer, indent=4)  # 增加缩进，不然可能得到的数据都在一行

print(answer_json)  # 打印查看

{
“distance”: 149,
“duration”: 127,
“direction”: 3,
“instruction”: “\u8d70150\u7c73,\u5230\u8fbe\u7ec8\u70b9”,
“path”: “116.46359089038,39.94356399372;116.4645233315,39.943548712251;116.46479282315,39.943537026419;116.46489702659,39.943532877602;116.46500302664,39.943527968169;116.46535067088,39.943512824986”,
“start_location”: {
“lng”: “116.46359089038”,
“lat”: “39.94356399372”
},
“end_location”: {
“lng”: “116.46535067088”,
“lat”: “39.943512824986”
}
}

通过以上的输出，可以看到，这里获取打印到的信息包括walking(步行导航)、riding(骑行导航)、driving(驾车导航)、transit(公交导航)以及
origin(起始点的经纬度坐标，前者为纬度，后者为经度)、destination(终点的经纬度坐标)等参数。
后续可根据相应的距离需求来进行获取。

瓦片块图像处理

针对地图图像处理，这里主要依赖于地图瓦片，简单理解就是根据不同地图层级的大小以及经纬度信息将地图图像根据行列划分为多个子图像，这个子图像就是地图瓦片。

瓦片地图模型概述及原理

瓦片地图金字塔模型是一种多分辨率层次模型，从瓦片金字塔的底层到顶层，分辨率越来越低，但表示的地理范围不变。首先确定地图服务平台所要提供的缩放级别的数量N，把缩放级别最高、地图比例尺最大的地图图片作为金字塔的底层，即第0层，并对其进行分块，从地图图片的左上角开始，从左至右、从上到下进行切割，分割成相同大小(比如256x256像素)的正方形地图瓦片，形成第0层瓦片矩阵;在第0层地图图片的基础上，按每像素分割为2×2个像素的方法生成第1层地图图片，并对其进行分块，分割成与下一层相同大小的正方形地图瓦片，形成第1层瓦片矩阵;采用同样的方法生成第2层瓦片矩阵;…;如此下去，直到第N一1层，构成整个瓦片金字塔。

其实切片之后的地图瓦片是栅格图像，并不具备定位信息，不过切片运用了相关切片算法之后，可以计算出具体定位的位置。例如采用WGS84大地坐标系为空间参考，对地图进行切片，采用一定的切片算法，例如用经纬度步长等比例分割形成地图瓦片，当需要对一个具体地方进行定位时，可以根据经纬度步长来计算具体位置，以此来达到定位的功能。

地图瓦片的基本原理

1.瓦片尺寸通常为256*256像素。

2.地图的最小zoom level是0，这时候，整个世界地图是一张瓦片，根据Web-Mercator投影公式可知，在zoom level=0的时候，经度[-180°,180°]，纬度[-85.051129°, 85.051129°]区间的数据可以被投影在一个256*256像素的正方形图片上。

3.zoom level越大，组成世界地图的瓦片数越多，可以展示的地图内容越详细。

4.某一zoom level下的地图瓦片是由它的上一层级的各瓦片切割成的4个瓦片组成，形成了瓦片金字塔。

5.瓦片通常都是png格式的，命名格式通常为z(zoom level)/x/y.png。

经维度坐标与瓦片的转换

经纬度坐标与瓦片的转换，最关键的参数是zoom level（z）。

瓦片z=0，x=0，y=0，对应的坐标范围是经度[-180,180]，纬度[-85.05,85.05]，对应的瓦片尺寸是256像素*246像素。

经纬度坐标(lng,lat)转瓦片坐标(tileX,tileY)的公式如下：

$tileX=int((lng+180)/360*2^z)$

$tileY=(1- arsinh(tan(lat*π/180))/π) *2^(z-1)$

其中z是zoom level，计算结果取整数。

其中反双曲正弦函数 $a rs inh (x)$ ，等价于 $ln(x+(x^2+1)^{0.5})$ 。

瓦片坐标(tileX,tileY)转经纬度坐标(lng,lat)公式，计算结果为该瓦片对应的最小经度和最大纬度：

$lng=tileX/2^z * 360-180$

$lat=arctan(sinh(π*(1-2*tileY/2^z)))*180/π$

其中，z是zoom level，双曲函数 $sinh(x)=(e^x-e^(-x))*{0.5}$ 。

经纬度坐标(lng,lat)转像素坐标(pixelX,pixelY)公式：

$pixelX=(lng+180)/360*2^z*256%256$

$pixelY=(1-ln(tan(lat*π/180)+sec(lat*π/180))/(2*π))*2^z*256%256$

经纬度、瓦片编号、像素坐标之间转换的python代码为：

import math
# 像素分辨率
def getResolution(level):
    return 156543.03*math.pow(2,-level)
# 经纬度转瓦片
def lnglatToTile(lng,lat,level):
    tileX=int((lng+180)/360*math.pow(2,level))
    tileY=int((1-math.asinh(math.tan(math.radians(lat)))/math.pi)*math.pow(2,level-1))
    return tileX,tileY
# 瓦片转经纬度
def tileToLnglat(tileX,tileY,level):
    lng=tileX/math.pow(2,level)*360-180
    lat=math.degrees(math.atan(math.sinh(math.pi*(1-2*tileY/math.pow(2,level)))))
    return lng,lat
# 经纬度转像素
def lnglatToPixel(lng,lat,level):
    pixelX=round((lng+180)/360*math.pow(2,level)*256%256)
    pixelY=round((1-math.log(math.tan(math.radians(lat))+1/math.cos(math.radians(lat)))/(2*math.pi))*math.pow(2,level)*256%256)
    return pixelX,pixelY
# 瓦片和像素转经纬度
def pixelToLnglat(tileX,tileY,pixelX,pixelY,level):
    lng=(tileX+pixelX/256)/math.pow(2,level)*360-180
    lat=math.degrees(math.atan(math.sinh(math.pi-2*math.pi*(tileY+pixelY/256)/math.pow(2,level))))
    return lng,lat
 
if __name__ == '__main__':
    # 北京坐标[116.391305,39.905530]
    print(getResolution(3))
    print(lnglatToTile(116.37,39.640,10))
    print(tileToLnglat(843, 388,10))
    print(lnglatToPixel(116.37,39.640,10))
    print(pixelToLnglat(843, 388,2,256,10))

19567.87875
(843, 388)
(116.3671875, 39.90973623453718)
(2, 256)
(116.36993408203125, 39.6395375643667)

根据经纬度批量获取下载瓦片图像

根据上面的操作与内容，已经知道了瓦片的基本概念含义以及如何去转换获取瓦片，
那么下面就基于百度地图API来进行相应经纬度瓦片图像的获取吧。
另外这里还需要说明的是，由于这里的地图简单分为两种，一种是抽象地图，一种是卫星地图，相应的瓦片也有不同的区分。
可以通过设置参数进行区分。详情可见注释。

这里先简单验证下，AK是可用的，也就是能够成功获取到信息。

import math
import urllib
import urllib.request as urllib2
import json
from urllib.request import urlopen
from urllib import parse

akey = 'dYg9wxs2v2xfaNEpTmcR4S6jZ0idkHN6'


def latlon2px(z, lat, lon):
    x = 2 ** z * (lon + 180) / 360 * 256
    y = -(.5 * math.log(
        (1 + math.sin(math.radians(lat))) / (1 - math.sin(math.radians(lat)))) / math.pi - 1) * 256 * 2 ** (z - 1)
    return x, y


def latlon2xy(z, lat, lon):
    x, y = latlon2px(z, lat, lon)
    x = int(x / 256)  # ,int(x%256)
    y = int(y / 256)  # ,int(y%256)
    return x, y


def bd_latlng2xy(z, lat, lng):
    # coords = str(lng) + ',' + str(lat)
    # url = 'http://api.map.baidu.com/geoconv/v1/?coords='+coords+'&ak=1jyKeVBnklxsB2Zyduy6wXdnwWEPUjaZ'
    url = 'http://api.map.baidu.com/geoconv/v1/?'
    # url = 'https://api.map.baidu.com/staticimage/v2/?'
    # url = 'https://api.map.baidu.com/staticimage/v2?ak=dYg9wxs2v2xfaNEpTmcR4S6jZ0idkHN6'
    args = {'coords': str(lng) + ',' + str(lat),
            'from': 5,
            'to': 6,
            'output': 'json',
            'ak': akey}
    data = urllib.parse.urlencode(args)
    response = urllib2.urlopen(url + data)
    print('response:', response)
    result = response.read()
    print('ress:', result)
    result = json.loads(result, encoding='utf-8')
    print('result:', result)
    loc = result["result"][0]
    res = 2 ** (18 - z)
    x = loc[u'x'] / res
    y = loc[u'y'] / res
    return x, y

z = 19
lat = 31.025819
lng = 121.434229
x, y = bd_latlng2xy(z, lat, lng)
print(x // 256)
print(y // 256)  # only right when lat>0 lng>0

response: <http.client.HTTPResponse object at 0x7f9a5c0faf50>
ress: b'{"status":0,"result":[{"x":13518143.623854652,"y":3613877.0373780529}]}'
result: {'status': 0, 'result': [{'x': 13518143.623854652, 'y': 3613877.0373780527}]}
105610.0
28233.0

批量下载瓦片图像-程序

import urllib.request as urllib2
import os, sys
import math

import time
import random


def download_tiles(zoom, lat_start, lat_stop, lon_start, lon_stop, satellite=True):
    start_x, start_y = bd_latlng2xy(zoom, lat_start, lon_start)
    stop_x, stop_y = bd_latlng2xy(zoom, lat_stop, lon_stop)

    start_x = int(start_x // 256)
    start_y = int(start_y // 256)
    stop_x = int(stop_x // 256)
    stop_y = int(stop_y // 256)

    print("x range", start_x, stop_x)
    print("y range", start_y, stop_y)
    root_save = "tiles"
    os.makedirs(root_save, exist_ok=True)

    for x in range(start_x, stop_x):
        for y in range(start_y, stop_y):

            url = None
            filename = None

            if satellite:
                url = "http://shangetu0.map.bdimg.com/it/u=x=%d;y=%d;z=%d;v=009;type=sate&fm=46&udt=20150504&app=webearth2&v=009&udt=20150601" % (
                x, y, zoom)
                filename = "%d_%d_%d_s.jpg" % (zoom, x, y)
            else:
                url = 'http://online3.map.bdimg.com/tile/?qt=tile&x=%d&y=%d&z=%d&styles=pl&scaler=1&udt=20180810' % (x, y, zoom)
                filename = "%d_%d_%d_r.png" % (zoom, x, y)
            filename = os.path.join(root_save, filename)
            print('filename', filename)
            if not os.path.exists(filename):

                bytes = None

                try:
                    req = urllib2.Request(url, data=None)
                    response = urllib2.urlopen(req)
                    bytes = response.read()
                except Exception as e:
                    print("--", filename, "->", e)
                    sys.exit(1)
                print("-- saving", filename)

                f = open(filename, 'wb')
                f.write(bytes)
                f.close()

                time.sleep(1 + random.random())

获取卫星地图瓦片

zoom = 19
lat_start, lon_start = 31.022547, 121.429391
lat_stop, lon_stop = 31.041453, 121.45749
satellite = True  # roads if false
download_tiles(zoom, lat_start, lat_stop, lon_start, lon_stop, satellite)

下面可以简单预览下。

import ipywidgets as widgets
import IPython.display as display
## Read images from file (because this is binary, maybe you can find how to use ByteIO) but this is more easy
img1 = open('/home/aistudio/tiles/19_105606_28236_s.jpg', 'rb').read()
img2 = open('/home/aistudio/tiles/19_105606_28238_s.jpg', 'rb').read()
img3 = open('/home/aistudio/tiles/19_105621_28238_s.jpg', 'rb').read()
## Create image widgets. You can use layout of ipywidgets only with widgets.
## Set image variable, image format and dimension.
wi1 = widgets.Image(value=img1, format='png', width=256, height=256)
wi2 = widgets.Image(value=img2, format='png', width=256, height=256)
wi3 = widgets.Image(value=img3, format='png', width=256, height=256)
## Side by side thanks to HBox widgets
sidebyside = widgets.HBox([wi1, wi2, wi3])
## Finally, show.
display.display(sidebyside)

HBox(children=(Image(value=b'\xff\xd8\xff\xe0\x00\x10JFIF\x00\x01\x01\x00\x00\x01\x00\x01\x00\x00\xff\xfe\x00>…

获取抽象地图瓦片

zoom = 19
lat_start, lon_start = 31.022547, 121.429391
lat_stop, lon_stop = 31.041453, 121.45749
satellite = False  # roads if false
download_tiles(zoom, lat_start, lat_stop, lon_start, lon_stop, satellite)

import ipywidgets as widgets
import IPython.display as display
## Read images from file (because this is binary, maybe you can find how to use ByteIO) but this is more easy
img1 = open('/home/aistudio/tiles/19_105606_28231_r.png', 'rb').read()
img2 = open('/home/aistudio/tiles/19_105606_28247_r.png', 'rb').read()
img3 = open('/home/aistudio/tiles/19_105606_28248_r.png', 'rb').read()
## Create image widgets. You can use layout of ipywidgets only with widgets.
## Set image variable, image format and dimension.
wi1 = widgets.Image(value=img1, format='png', width=256, height=256)
wi2 = widgets.Image(value=img2, format='png', width=256, height=256)
wi3 = widgets.Image(value=img3, format='png', width=256, height=256)
## Side by side thanks to HBox widgets
sidebyside = widgets.HBox([wi1, wi2, wi3])
## Finally, show.
display.display(sidebyside)

HBox(children=(Image(value=b'\x89PNG\r\n\x1a\n\x00\x00\x00\rIHDR\x00\x00\x01\x00\x00\x00\x01\x00\x08\x03\x00\x…

合并瓦片地图

上面已经分别获取到了卫星地图和抽象地图在指定经纬度的256*256的瓦片地图，那么怎么样得到完整区域下的地图呢，
可以通过下面程序将其进行行列瓦片图的合并。

from PIL import Image
import sys, os

def merge_tiles(zoom, lat_start, lat_stop, lon_start, lon_stop, satellite=True):
    TYPE, ext = 'r', 'png'
    if satellite:
        TYPE, ext = 's', 'jpg'

    x_start, y_start = bd_latlng2xy(zoom, lat_start, lon_start)
    x_stop, y_stop = bd_latlng2xy(zoom, lat_stop, lon_stop)
    x_start = int(x_start / 256)
    y_start = int(y_start / 256)
    x_stop = int(x_stop / 256)
    y_stop = int(y_stop / 256)
    print("x range", x_start, x_stop)
    print("y range", y_start, y_stop)

    w = (x_stop - x_start) * 256
    h = (y_stop - y_start) * 256

    print("width:", w)
    print("height:", h)

    result = Image.new("RGBA", (w, h))
    root_save = "tiles"

    print("start download===",x_start, x_stop,y_start,y_stop)
    for x in range(x_start, x_stop):
        for y in range(y_start, y_stop):

            filename = "%d_%d_%d_%s.%s" % (zoom, x, y, TYPE, ext)
            filename = os.path.join(root_save, filename)

            if not os.path.exists(filename):
                print("-- missing", filename)
                continue

            x_paste = (x - x_start) * 256
            y_paste = (y_stop - y - 1) * 256

            try:
                i = Image.open(filename)
            except Exception as e:
                print("-- %s, removing %s" % (e, filename))
                trash_dst = os.path.expanduser("~/.Trash/%s" % filename)
                os.rename(filename, trash_dst)
                continue

            result.paste(i, (x_paste, y_paste))

            del i
    result = result.convert("RGB")
    result.save("map_%s.%s" % (TYPE, ext))

合并卫星地图瓦片

zoom = 19
lat_start, lon_start = 31.022547, 121.429391
lat_stop, lon_stop = 31.041453, 121.45749
satellite = True  # roads if false
# zoom = 15
# lat_start, lon_start = 29.373026, 106.25131
# lat_stop, lon_stop = 29.755033, 106.72734
merge_tiles(zoom, lat_start, lat_stop, lon_start, lon_stop, satellite=True)

response: <http.client.HTTPResponse object at 0x7f0db864da50>
ress: b'{"status":0,"result":[{"x":13517605.054298275,"y":3613454.1099721884}]}'
result: {'status': 0, 'result': [{'x': 13517605.054298274, 'y': 3613454.1099721882}]}
response: <http.client.HTTPResponse object at 0x7f0d99269cd0>
ress: b'{"status":0,"result":[{"x":13520733.05470436,"y":3615898.0367417645}]}'
result: {'status': 0, 'result': [{'x': 13520733.05470436, 'y': 3615898.0367417643}]}
x range 105606 105630
y range 28230 28249
width: 6144
height: 4864
start download=== 105606 105630 28230 28249

合并完成后，进行效果预览下。

%matplotlib inline
def imshow_img(img_path):

    import matplotlib.pyplot as plt
    import matplotlib.image as mpimg
    
    img = mpimg.imread(img_path)
    plt.figure(figsize=(10,10))
    plt.imshow(img)
    plt.axis('off')
    plt.show()
img_path = '/home/aistudio/map_s.jpg'
imshow_img(img_path)

在这里插入图片描述

合并抽象地图瓦片

zoom = 19
lat_start, lon_start = 31.022547, 121.429391
lat_stop, lon_stop = 31.041453, 121.45749
satellite = True  # roads if false
merge_tiles(zoom, lat_start, lat_stop, lon_start, lon_stop, satellite=False)

response: <http.client.HTTPResponse object at 0x7f9a458032d0>
ress: b'{"status":0,"result":[{"x":13517605.054298275,"y":3613454.1099721884}]}'
result: {'status': 0, 'result': [{'x': 13517605.054298274, 'y': 3613454.1099721882}]}
response: <http.client.HTTPResponse object at 0x7f9a5c0bc050>
ress: b'{"status":0,"result":[{"x":13520733.05470436,"y":3615898.0367417645}]}'
result: {'status': 0, 'result': [{'x': 13520733.05470436, 'y': 3615898.0367417643}]}
x range 105606 105630
y range 28230 28249
width: 6144
height: 4864
start download=== 105606 105630 28230 28249

ad=== 105606 105630 28230 28249

img_path = '/home/aistudio/map_r.png'
imshow_img(img_path)

图像分类处理

(!!!下面代码需要在本地环境运行，目前AI Studio环境暂不支持)

处理完以上步骤后，可以将图片从AI Studio上打包下载下来，然后运行下面的代码，则可以通过按键设置来进行图片预览和分类。
这里简单的将卫星瓦片图像分为房屋、道路以及空地三个类。
效果如下，后续可通过这种方式来制作卫星瓦片地图的分类数据集，同时下面的手动图像分类代码也可用于其他图像的分类。

import cv2
import os
import shutil
import time
import sys

leftkeys = (81, 110, 65361, 2424832)
rightkeys = (83, 109, 65363, 2555904)
# 当前脚本工作的目录路径
root_dir = os.getcwd()
# print(root_dir)
# os.path.abspath()获得绝对路径
root_absdir = os.path.abspath(os.path.dirname(__file__))
print(root_absdir)

def make_dirs2(n):
    if not os.path.exists(os.path.join(root_dir, n)):
        os.makedirs(os.path.join(root_dir, n))


def Classification_Tools(data_dir,class_names,action_type='copy'):

    global image
    if not os.path.exists(data_dir):
        print('data_all not exists, please put data to: ', data_dir)
        time.sleep(5)
        exit()

    unconfirmed = 'unconfirmed'  # 不确定数据存放路径
    make_dirs2(unconfirmed)


    num_cls = len(class_names)
    for j in range(num_cls):
        make_dirs2(str(class_names[j]))

    image_list = os.listdir(data_dir)
    if len(image_list) == 0:
        print('no image in %s ... please put data to: %s' % (data_dir, data_dir))
        time.sleep(5)
        exit()
    cv2.namedWindow('Classification_Tools', 0)
    i = 0
    coccus_label = None
    for i in range(0, len(image_list)):
        assert i < len(image_list), ('no image left...')
        print('i', i)
        image_path = os.path.join(data_dir, image_list[i])
        print(image_path)
        image = cv2.imread(image_path)
        if image.shape is None:
            print('The images open fail!')
            continue
        else:
            print(image.shape)
            cv2.imshow('Classification_Tools', image)
            key = cv2.waitKeyEx()

        if key == ord('d'):
            coccus_label = 'unconfirmed'
            if action_type == 'copy':
                shutil.copy(image_path, os.path.join(root_dir, unconfirmed))
            else:
                shutil.move(image_path, os.path.join(root_dir, unconfirmed))
            i += 1  # (i + 1) % len(image_list)
        if key in rightkeys:
            i += 1  # (i + 1) % len(image_list)
        if key in leftkeys and coccus_label != None:
            print('leftkeys:', os.path.join(('./' + str(coccus_label)), image_list[i - 1]))
            if os.path.exists(os.path.join(('./' + str(coccus_label)), image_list[i - 1])):
                print('successful', os.path.join(('./' + str(coccus_label)), image_list[i - 1]))
                if action_type == 'copy':
                    shutil.copy(os.path.join(('./' + str(coccus_label)), image_list[i - 1]), data_dir)
                else:
                    shutil.move(os.path.join(('./' + str(coccus_label)), image_list[i - 1]), data_dir)
            i -= 1
            if i < 0:
                i = len(image_list) - 1


        if (key == ord('q')) or (key == 27):
            break

        for j in range(num_cls):
            if key & 0xFF == ord(str(j)):
                coccus_label = str(j)
                shutil.move(image_path, os.path.join(root_dir, str(class_names[j])))

                i += 1  # (i + 1) % len(image_list)
                break

if __name__ == '__main__':
    # 待分类数据路径
    data_dir = '/home/linxu/Projects/PycharmProjects/AI-Map/core/tiles'
    # 待分类名称列表
    class_names = ['house', 'road', 'space']
    # 分类标注器
    Classification_Tools(data_dir,class_names,)