python爬虫,爬取遥感影像瓦片并自动拼接

效果

实现思路

首先根据经纬度范围确定能够下载的瓦片最高等级(影像最大分辨率),然后下载瓦片,将瓦片按照顺序拼接起来 

 源代码

urllib3 == 1.26.7

requests == 2.26.0

numpy == 1.21.4

Pillow == 8.4.0

opencv-python == 4.5.4.58

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 29 16:41:43 2021

@author: yangzhen
2021.4.29
"""


from urllib import request
import re
import urllib.request
import random
import math
import os
import numpy as np
import io
import PIL.Image as pil
import cv2


def deg2num(lat_deg, lon_deg, zoom):
    """经纬度反算切片行列号 3857坐标系"""
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)


def getimg(Tpath, x, y):
    """下载瓦片"""
    # 简单的反爬虫
    agents = [
    'Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/60.0.3112.101 Safari/537.36',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/532.5 (KHTML, like Gecko) Chrome/4.0.249.0 Safari/532.5',
    'Mozilla/5.0 (Windows; U; Windows NT 5.2; en-US) AppleWebKit/532.9 (KHTML, like Gecko) Chrome/5.0.310.0 Safari/532.9',
    'Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US) AppleWebKit/534.7 (KHTML, like Gecko) Chrome/7.0.514.0 Safari/534.7',
    'Mozilla/5.0 (Windows; U; Windows NT 6.0; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/9.0.601.0 Safari/534.14',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/10.0.601.0 Safari/534.14',
    'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.20 (KHTML, like Gecko) Chrome/11.0.672.2 Safari/534.20", "Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/534.27 (KHTML, like Gecko) Chrome/12.0.712.0 Safari/534.27',
    'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.1 (KHTML, like Gecko) Chrome/13.0.782.24 Safari/535.1']
    try:
        print(str(x) + '_' + str(y) + '下载成功')
        req = urllib.request.Request(Tpath)
        req.add_header('User-Agent', random.choice(agents))  # 换用随机的请求头
        pic = urllib.request.urlopen(req, timeout=60)
        timg = pic.read()
        outimg = pil.new('RGBA', (256, 256))
        outimgfile = io.BytesIO(timg)
        img = pil.open(outimgfile)
        img = cv2.cvtColor(np.asarray(img), cv2.COLOR_RGB2BGR)
    except Exception:
        print(str(x) + '_' + str(y) + '下载失败,重试')
        getimg(Tpath, x, y)
    return img
        
def GetMapTile(lat1, lon1, lat2, lon2):
    """根据经纬度范围获取背景"""
    zooms = []
    for i in range(1, 19):
        l = deg2num(lat1, lon1, i)
        r = deg2num(lat2, lon2, i)
        if l[0] - r[0] == 0 or l[1] - r[1] == 0:
            continue
        else:
            zooms.append(i)
    # 根据有瓦片数据的地方下载瓦片
    zoom = zooms[-1]
    # 根据经纬度确定瓦片位置
    lefttop = deg2num(lat1, lon1, zoom)
    rightbottom = deg2num(lat2, lon2, zoom)
    # 一边下载一边拼接
    imgcolumns = []
    imglist_all = []
    for x in range(lefttop[0], rightbottom[0]):
        imgrows = []
        imglist_row = []
        for y in range(lefttop[1], rightbottom[1]):
            #Google地图瓦片
            # tilepath = 'http://www.google.cn/maps/vt/pb=!1m4!1m3!1i'+str(zoom)+'!2i'+str(x)+'!3i'+str(y)+'!2m3!1e0!2sm!3i345013117!3m8!2szh-CN!3scn!5e1105!12m4!1e68!2m2!1sset!2sRoadmap!4e0'
            #Google影像瓦片
            # tilepath = 'http://mt3.google.cn/vt/lyrs=s@110&hl=zh-CN&gl=cn&src=app&x='+str(x)+'&y='+str(y)+'&z='+str(zoom)+'&s=G'
            #天地图-地图
            #tilepath = 'http://t4.tianditu.com/DataServer?T=vec_w&x='+str(x)+'&y='+str(y)+'&l='+str(zoom)+'&tk=45c78b2bc2ecfa2b35a3e4e454ada5ce'
            #天地图-标注
            #tilepath = 'http://t3.tianditu.com/DataServer?T=cva_w&x='+str(x)+'&y='+str(y)+'&l='+str(zoom)+'&tk=45c78b2bc2ecfa2b35a3e4e454ada5ce'
            #天地图-影像
            # tilepath = 'http://t2.tianditu.gov.cn/DataServer?T=img_w&x='+str(x)+'&y='+str(y)+'&l='+str(zoom)+'&tk=2ce94f67e58faa24beb7cb8a09780552'
            #天地图-影像标注
            # tilepath = 'http://t2.tianditu.gov.cn/DataServer?T=cia_w&x='+str(x)+'&y='+str(y)+'&l='+str(zoom)+'&tk=2ce94f67e58faa24beb7cb8a09780552'
            #高德地图影像瓦片
            tilepath = "http://wprd01.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scl=1&style=6&x=" + \
                        str(x) + "&y=" + str(y) + "&z=" + str(zoom) + "&ltype=3"
            img = getimg(tilepath, x, y)
            imgrows.append(img)
        imgcolumn = np.vstack(imgrows)
        imgcolumns.append(imgcolumn)
    finalimg = np.hstack(imgcolumns)
    return finalimg


if __name__ == '__main__':
    img = GetMapTile(28.173161534026335,112.92222966168211,28.164762736523976,112.9360054871826)
    cv2.namedWindow('test', 0)
    cv2.imshow('test', img)
    cv2.waitKey()

 

  • 3
    点赞
  • 14
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值