谷歌卫星图瓦片介绍

1
2
http://mt{index}.google.com/vt/lyrs={style}&hl=zh-CN&gl=CN&src=app&x={x}&y={y}&z={z}"

参数解析

  • {index}为 0 到 3，分别为谷歌的 4 个服务器，随便取一个就行~

  • {style}为地图的类型，比如卫星图、路线图什么的.

    • m：路线图
    • t：地形图
    • p：带标签的地形图
    • s：卫星图
    • y：带标签的卫星图
    • h：标签层（路名、地名等）
    • 其中 m 路线图和 s 卫星图是我们想要的。

  • {x}{y}{z} 则分别是瓦片坐标的 xy 和缩放级别 z，z 取 0-22。但是我测试发现卫星图最大只能取到 20。不过即使是路线图，到 20 级也就足够用了。

  • gl 参数决定了卫星图是否会偏移

    • 假设不加这个参数，会按 WGS-84 来，会和路线图有偏移（因为路线图只有 GCJ-02 的版本）
    • 加了这个参数（gl=CN）就会使卫星图也变成 GCJ-02。
    • 也就是说，如果你使用 www.google.cn/maps 或者 ditu.google.cn 来获取坐标，那么最后你还按照这个坐标进行瓦片的获取，同时加上 gl 参数即可。如果你想按照 WGS-84 来获取坐标，那么就访问国际版的谷歌地图，获取瓦片时不加 gl 参数就好了。

对于 python3，可以用 urllib.request 库进行图片的下载，然后用 PIL.Image 库（pillow）进行图片的合并。

下载关键代码

伪墨卡托与经纬度互转

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#------------------wgs84与web墨卡托互转-------------------------

# WGS-84经纬度转Web墨卡托
def wgs_to_macator(x, y):
    y = 85.0511287798 if y > 85.0511287798 else y
    y = -85.0511287798 if y < -85.0511287798 else y

    x2 = x * 20037508.34 / 180
    y2 = log(tan((90 + y) * pi / 360)) / (pi / 180)
    y2 = y2 * 20037508.34 / 180
    return x2, y2

# Web墨卡托转经纬度
def mecator_to_wgs(x, y):
    x2 = x / 20037508.34 * 180
    y2 = y / 20037508.34 * 180
    y2 = 180 / pi * (2 * atan(exp(y2 * pi / 180)) - pi / 2)
    return x2, y2

#-------------------------------------------------------------

根据经纬度计算瓦片位置

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
# 根据WGS-84 的经纬度获取谷歌地图中的瓦片坐标
def wgs84_to_tile(j, w, z):
    '''
    Get google-style tile cooridinate from geographical coordinate
    j : Longittude
    w : Latitude
    z : zoom
    '''
    isnum = lambda x: isinstance(x, int) or isinstance(x, float)
    if not(isnum(j) and isnum(w)):
        raise TypeError("j and w must be int or float!")

    if not isinstance(z, int) or z < 0 or z > 22:
        raise TypeError("z must be int and between 0 to 22.")

    if j < 0:
        j = 180 + j
    else:
        j += 180
    j /= 360  # make j to (0,1)

    w = 85.0511287798 if w > 85.0511287798 else w
    w = -85.0511287798 if w < -85.0511287798 else w
    w = log(tan((90 + w) * pi / 360)) / (pi / 180)
    w /= 180  # make w to (-1,1)
    w = 1 - (w + 1) / 2  # make w to (0,1) and left top is 0-point

    num = 2**z
    x = floor(j * num)
    y = floor(w * num)
    return x, y


经纬度度分秒与小数互转

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
def dms2dd(degrees, minutes, seconds, direction):

    dd = float(degrees) + float(minutes)/60 + float(seconds)/(60*60);

    if direction == 'E' or direction == 'N':

        dd *= -1

    return dd;

def dd2dms(deg):

    d = int(deg)

    md = abs(deg - d) * 60

    m = int(md)

    sd = (md - m) * 60

    return [d, m, sd]

def parse_dms(dms):

    parts = re.split('[^\d\w]+', dms)

    lat = dms2dd(parts[0], parts[1], parts[2], parts[3])

    return (lat)

下载器

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
class Downloader(Thread):
    # multiple threads downloader
    def __init__(self,index,count,urls,datas,update):
        # index 表示第几个线程，count 表示线程的总数，urls 代表需要下载url列表，datas代表要返回的数据列表。
        # update 表示每下载一个成功就进行的回调函数。
        super().__init__()
        self.urls=urls
        self.datas=datas
        self.index=index
        self.count=count
        self.update=update

    def download(self,url):
        HEADERS = {'User-Agent': random.choice(agents)}
        # print(url)

        err=0
        while(err<4):
            full_url = url.replace('_index_',str(err))  # 如果下载失败，则切换服务器
            header = ur.Request(full_url,headers=HEADERS)
            try:
                data = ur.urlopen(header).read()
            except:
                err+=1
            else:
                return data
        # raise Exception("Bad network link.")
        print("Bad network link.",url)

    def run(self):
        for i,url in enumerate(self.urls):
            if i%self.count != self.index:
                continue
            self.datas[i]=self.download(url)
            if mutex.acquire():
                self.update()
                mutex.release()

读取shapefile文件

读取shapefile文件的方式比较多，这里使用了两种方式，一种是使用shapefile模块读取，一种是使用gdal读取

shapefile读取

1
2
3
4
shp = shapefile.Reader(shapefilename)
border_shapes = shp.shapes()
area_infors = shp.records()
return border_shapes, area_infors

gdal读取

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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def getShpFileByGDAL(shapefile, query_name):
    '''
    结果保存方式为b,c,h,w
    b为有多少个区域
    c为区域轮廓个数，例如c=1，表示只有一个外轮廓，无空洞；c=2，表示下标为0的位置为外轮廓，下标为1的位置为外轮廓，除了0以外都是外轮廓
    h为轮廓点的个数
    w为2，表示经纬度
    '''
    # 支持中文路径
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
    # 支持中文编码
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 注册所有的驱动
    ogr.RegisterAll()
    # 打开数据
    print('打开', shapefile)
    ds = ogr.Open(shapefile, 0)
    if ds == None:
        return "打开文件失败！"

    # 获取数据源中的图层个数，shp数据图层只有一个，gdb、dxf会有多个
    iLayerCount = ds.GetLayerCount()
    print("\t图层个数 = ", iLayerCount)
    # 获取第一个图层
    result = []
    bbox = [10000, 10000, -10000, -10000]
    for layerIdx in range(iLayerCount):
        oLayer = ds.GetLayerByIndex(layerIdx)
        if oLayer == None:
            return "获取图层失败！"
        # 对图层进行初始化
        oLayer.ResetReading()
        # 输出图层中的要素个数
        num = oLayer.GetFeatureCount(0)
        print("\t要素个数 = ", num)
        result_list = []
        # 获取要素
        for i in range(0, num):
            ofeature = oLayer.GetFeature(i)
            # id = ofeature.GetFieldAsString("id")
            name = ofeature.GetFieldAsString('name')

            if query_name in name:

                count = ofeature.GetGeometryRef().GetGeometryCount()
                for gemIdx in range(count):
                    gemo = ofeature.GetGeometryRef().GetGeometryRef(gemIdx)
                    gemo_data = gemo.ExportToJson()
                    gemo_json = json.loads(gemo_data)
                    gemo_np = np.array(gemo_json['coordinates'])
                    print(name, gemo_np.shape)

                    min_x = min(np.min(gemo_np[..., 0]), bbox[0])
                    max_x = max(np.max(gemo_np[..., 0]), bbox[2])
                    min_y = min(np.min(gemo_np[..., 1]), bbox[1])
                    max_y = max(np.max(gemo_np[..., 1]), bbox[3])

                    bbox = [min_x, min_y, max_x, max_y]
                    result.append(gemo_np)
                # print(dir(ofeature.GetGeometryRef()),
                #       ofeature.GetGeometryRef().Boundary())

    ds.Destroy()
    del ds
    return result, bbox

裁剪

最后的裁剪使用的是用一张mask对图片做掩码裁剪

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 制作掩码
mask = np.zeros((height, width), np.uint8)

for points_list in mercator_list:
    num = len(points_list)
    for idx in range(num):
        fill_value = 0
        if idx == 0:
            fill_value = 1
        else:
            fill_value = 0
        mask = cv2.fillPoly(mask, [np.array(points_list[idx])], fill_value)

效果
image

完整代码

参考博客

https://blog.csdn.net/mygisforum/article/details/22997879
https://www.cnblogs.com/modyuan/p/12293843.html
https://blog.csdn.net/u012866178/article/details/83188581
https://www.cnblogs.com/jyughynj/p/pythonGetMapTile.html
https://www.cnblogs.com/xinwenpeng/p/9818953.html