【无标题】python自动下载瓦片影像

利用python下载瓦片影像，带坐标，适合下载小范围，大范围影像数据会被服务器阻挡

import math
import multiprocessing
import numpy as np
import os
import queue
import shutil
import threading
from PIL import Image
from datetime import datetime
from urllib import request
import time
import re
from osgeo import gdal,osr

class MapDownloader(object):
    def __init__(self, lat_start, lng_start, lat_end, lng_end, zoom=12, tile_size=512):

        self.tile_server = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}'
        self.lat_start = lat_start
        self.lng_start = lng_start
        self.lat_end = lat_end
        self.lng_end = lng_end
        self.zoom = zoom
        self.tile_size = tile_size

        self.q = queue.Queue()
        self.num_worker = multiprocessing.cpu_count()#线程数，如果线程过多，容易出现访问服务器被拒绝
        # print('线程数',self.num_worker)
        self._generate_xy_point()

    def _generate_xy_point(self):
        self._x_start, self._y_start = self._convert_latlon_to_xy(self.lat_start, self.lng_start)
        self._x_end, self._y_end = self._convert_latlon_to_xy(self.lat_end, self.lng_end)

    def _convert_latlon_to_xy(self, lat, lng):
        tiles_count = 1 << self.zoom

        point_x = (self.tile_size / 2 + lng * self.tile_size / 360.0) * tiles_count // self.tile_size
        sin_y = math.sin(lat * (math.pi / 180.0))
        point_y = ((self.tile_size / 2) + 0.5 * math.log((1 + sin_y) / (1 - sin_y)) *
                   -(self.tile_size / (2 * math.pi))) * tiles_count // self.tile_size

        return int(point_x), int(point_y)

    def _fetch_worker(self):
        while True:
            item = self.q.get()
            if item is None:
                break
            idx, url, current_tile = item
            # time.sleep(2)
            request.urlretrieve(url, current_tile)
            self.q.task_done()

    def write_into(self, filename):
        # 创建底图瓦片存储路径
        directory = os.path.abspath('./{}'.format(datetime.now().strftime("%Y-%m-%d_%H-%M-%S")))
        if not os.path.exists(directory):
            os.makedirs(directory)

        # generate source list
        idx = 1
        for x in range(0, self._x_end + 1 - self._x_start):
            for y in range(0, self._y_end + 1 - self._y_start):
                url = self.tile_server.format(
                    x=str(self._x_start + x), y=str(self._y_start + y), z=str(self.zoom))
                current_tile = os.path.join(directory, 'tile-{}_{}_{}.tif'.format(
                    str(self._x_start + x), str(self._y_start + y), str(self.zoom)))
                self.q.put((idx, url, current_tile))
                idx += 1

        # stop workers
        for i in range(self.num_worker):
            self.q.put(None)
        # 开始使用多线程来获取图块以加快处理
        self.q_size = self.q.qsize()
        threads = []
        for i in range(self.num_worker):
            t = threading.Thread(target=self._fetch_worker)
            t.start()
            threads.append(t)
        for t in threads:
            t.join()
        # 将图像合并为单个
        width, height = 256 * (self._x_end + 1 - self._x_start), 256 * (self._y_end + 1 - self._y_start)
        map_img = Image.new('RGB', (width, height))

        for x in range(0, self._x_end + 1 - self._x_start):
            for y in range(0, self._y_end + 1 - self._y_start):
                current_tile = os.path.join(directory, 'tile-{}_{}_{}.tif'.format(
                    str(self._x_start + x), str(self._y_start + y), str(self.zoom)))
                im = Image.open(current_tile)
                map_img.paste(im, (x * 256, y * 256))

        map_img.save(filename)

        # remove temp dir
        # shutil.rmtree(directory)
def main(lat_start, log_start, lat_end, log_end,zoomlevel,tile_size,imagename):
    try:
        md = MapDownloader(lat_start, log_start, lat_end, log_end,zoomlevel,tile_size)
        md.write_into(imagename+'.tif')
        print("已完成")
    except Exception as e:
        print("尝试创建image. Cause: {}".format(e))
def getSRSPair(dataset):
    '''
    获得给定数据的投影参考系和地理参考系
    :param dataset: GDAL地理数据
    :return: 投影参考系和地理参考系
    '''
    prosrs = osr.SpatialReference()
    prosrs.ImportFromWkt(dataset.GetProjection())
    geosrs = prosrs.CloneGeogCS()
    return prosrs, geosrs

def geo2lonlat(dataset, x, y):
    '''
    将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
    :param dataset: GDAL地理数据
    :param x: 投影坐标x
    :param y: 投影坐标y
    :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(prosrs, geosrs)
    coords = ct.TransformPoint(x, y)
    return coords[:2]

def def_geoCoordSys(read_path, img_transf, img_proj):
        array_dataset = gdal.Open(read_path)
        img_array = array_dataset.ReadAsArray(
            0, 0, array_dataset.RasterXSize, array_dataset.RasterYSize)
        if 'int8' in img_array.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in img_array.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        if len(img_array.shape) == 3:
            img_bands, im_height, im_width = img_array.shape
        else:
            img_bands, (im_height, im_width) = 1, img_array.shape

        filename = read_path[:-4] + '_proj' + '.tif'
        driver = gdal.GetDriverByName("GTiff")  # 创建文件驱动
        dataset = driver.Create(filename, im_width, im_height, img_bands, datatype)
        dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数
        dataset.SetProjection(img_proj)  # 写入投影
        # 写入影像数据
        if img_bands == 1:
            dataset.GetRasterBand(1).WriteArray(img_array)
        else:
            for i in range(img_bands):
                dataset.GetRasterBand(i + 1).WriteArray(img_array[i])
        print('文件路径',read_path, '坐标赋予成功!')

if __name__ == '__main__':
    starttime=time.time()
    regionpath='对应范围的影像路径（须有地理坐标）'
    zoomlevel=16 #影像等级
    imagename='Google_'
    dataset1 = gdal.Open(regionpath)
    if dataset1 is None:
        print('数据未读入')
        exit(0) 
    geom=dataset1.GetGeoTransform()
    row1=dataset1.RasterXSize     #行大小
    col1=dataset1.RasterYSize     #列大小
    x1=geom[0]
    y1=geom[3]
    log_start,lat_start=geo2lonlat(dataset1,geom[0],geom[3])  #左上角  返回经度、纬度

    x3=x1+row1*geom[1]
    y3=y1+col1*geom[5]
    log_end,lat_end= geo2lonlat(dataset1, x3, y3)             #右下角  返回经度、纬度
    print('范围经纬度',lat_start,log_start,lat_end,log_end)
    
    main(lat_start,log_start,lat_end,log_end,zoomlevel,512,imagename+str(zoomlevel))   #下载主函数输入参数  '纬度起点Latitude''经度起点Longitude''纬度终点''经度终点'

    image_no='E:/Vscode/'+imagename+str(zoomlevel)+'.tif'    #输出路径
    dataset2=gdal.Open(image_no)
    row2=dataset2.RasterXSize     #行大小
    col2=dataset2.RasterYSize     #列大小
    X_Resolution=(x3-x1)/row2
    Y_Resolution=(y3-y1)/col2
    img_pos_transf=(x1,X_Resolution , 0.0, y1, 0.0,Y_Resolution )
    print('Google坐标信息:',img_pos_transf)
    img_pos_proj=dataset1.GetProjection()
    def_geoCoordSys(image_no, img_pos_transf, img_pos_proj)   #坐标转换的函数
    print('下载所需时间',(time.time()-starttime))