【GIS】自动生成带坐标信息的拼接卫星影像图片（TIFF格式）

约定不等于承诺〃 2022-10-22 04:13 945阅读 0赞

### 文章目录 ###

*  0 项目需求
 *  1 实现思路
 *   *  1.1 技术流程
     *  1.2 用什么数据源？
     *  1.3 如何为影像增加坐标信息？
 *  2 具体实现
 *   *  2.1 下载卫星影像瓦片数据
     *  2.2 拼接卫星影像瓦片数据
     *  2.3 增加坐标信息
 *  3 代码封装
 *   *  3.1 感兴趣区域获取
     *  3.2 经纬度行列号转换
     *  3.3 Flask App

# 0 项目需求 #

经常会有用到卫星影像作为参考底图的时候，想要获取任一区域的卫星影像，可以借助`LSV`、`BigeMap`、奥维地图、`91`地图助手等地图软件，但这些软件有一些问题，要不然因为图源使用的是`Google Maps`而使得连接不稳定、数据经常不可下载，要不然就是对大比例尺卫星影像下载有限制，要求收费或者免费但是会打水印。

其实，前述软件自己并不生产卫星影像，他们做的工作不过是把图商（`Google Maps`、百度底图、高德底图、天地图）影像数据做了拼接并额外添加了地理信息而已。既然思路并不复杂，何不考虑自己实现之，毕竟“自己动手，丰衣足食”嘛。

# 1 实现思路 #

## 1.1 技术流程 ##

要把地图商业公司的瓦片数据变成带坐标信息的影像需要几步？答案是3步：

1.  **下载卫星影像瓦片数据**：图商瓦片一般是`256像素x256像素`大小的`JPG`或`PNG`格式的图片，按照缩放级别（`ZOOM`），行号（`ROW`）、列号（`COL`）形式组织（不同图商组织逻辑有一定出入，比如参考原点、级别数量），浏览某个区域的时候其实是用缩放级别、行号、列号组装图片`URI`，下载图片后在前端拼装到指定位置即实现我们查看地图时候效果。可以参考瓦片组织形式，组装瓦片图片url，请求后下载到本地；
2.  **拼接卫星影像瓦片数据**：我们需要的不是零散的瓦片，而是拼接好的带有坐标信息的影像图片，因此，必须对下载好后的所有瓦片数据做拼接处理；
3.  **增加坐标信息**：很显然，拼接好的`JPG`或者`PNG`无法保留坐标信息，而地理信息领域中`TIFF`格式是保存地理信息的常见格式，可以从瓦片的行列号信息中反算坐标信息，并将坐标信息保存至TIFF格式中。

## 1.2 用什么数据源？ ##

*  `Google Maps`的影像数据清晰、缩放级别多，是理想数据源，但因为网络问题无法访问（除非自己搭梯子），不能作为数据源的选项。
 *  百度地图影像数据尚可，但百度地图采用的是自己二次加偏的`BD09`坐标系，与`WGS84`或者`CGCS2000`相比，空间上有较大偏差，非得使用的话还得增加额外的投影变换步骤，略麻烦。
 *  天地图影像数据采用的是`CGCS2000`坐标系，有`1-18`个缩放级别，影像较为清晰，凑合能用，选择。

## 1.3 如何为影像增加坐标信息？ ##

`GDAL`帮我们实现了这个功能。只要有一张图片和这张图片四角坐标信息，`GDAL`就能将其转为带坐标信息的TIFF格式图片。

> GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。OGR是GDAL项目的一个分支，功能与GDAL类似，只不过它提供对矢量数据的支持。有很多著名的GIS类产品都使用了GDAL/OGR库，包括ESRI的ARCGIS 9.3，Google Earth和跨平台的GRASS GIS系统。

# 2 具体实现 #

## 2.1 下载卫星影像瓦片数据 ##

给定行列号范围，组装出瓦片URI，请求后将瓦片图片按照`ROW_COL.jpg`格式保存至同一个文件夹中。基于`Python urllib`库实现此功能。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JhaWR1XzI2NjQ2MTI5_size_16_color_FFFFFF_t_70]

示例代码如下：

# -*- coding: utf-8 -*-
    import http, time, random
    import urllib.request
    
    def getHeader():
    	""" 随机请求头 """
        USER_AGENTS = [
            "Mozilla/5.0 (compatible; MSIE 9.0; Windows NT 6.1; Win64; x64; Trident/5.0; .NET CLR 3.5.30729; .NET CLR 3.0.30729; .NET CLR 2.0.50727; Media Center PC 6.0)",
            "Mozilla/5.0 (compatible; MSIE 8.0; Windows NT 6.0; Trident/4.0; WOW64; Trident/4.0; SLCC2; .NET CLR 2.0.50727; .NET CLR 3.5.30729; .NET CLR 3.0.30729; .NET CLR 1.0.3705; .NET CLR 1.1.4322)",
            "Mozilla/4.0 (compatible; MSIE 7.0b; Windows NT 5.2; .NET CLR 1.1.4322; .NET CLR 2.0.50727; InfoPath.2; .NET CLR 3.0.04506.30)",
            "Mozilla/5.0 (Windows; U; Windows NT 5.1; zh-CN) AppleWebKit/523.15 (KHTML, like Gecko, Safari/419.3) Arora/0.3 (Change: 287 c9dfb30)",
            "Mozilla/5.0 (X11; U; Linux; en-US) AppleWebKit/527+ (KHTML, like Gecko, Safari/419.3) Arora/0.6",
            "Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv:1.8.1.2pre) Gecko/20070215 K-Ninja/2.1.1",
            "Mozilla/5.0 (Windows; U; Windows NT 5.1; zh-CN; rv:1.9) Gecko/20080705 Firefox/3.0 Kapiko/3.0",
            "Mozilla/5.0 (X11; Linux i686; U;) Gecko/20070322 Kazehakase/0.4.5"
        ]
    
        user_agent = random.choice(USER_AGENTS)
        headers = { 
            'User-Agent':user_agent
        }
    
        return headers
       
     
    def requestImg(url, file_name, num_retries=3):
    	""" 下载并保存图片 url：图片URI file_name：保存文件名 num_retries：失败重连次数 """
        try:
            headers = getHeader()
            request = urllib.request.Request(url, headers=headers)
            res = urllib.request.urlopen(request)
            with open(file_name, "wb") as f:
                content = res.read()
                f.write(content)
                res.close()
        except urllib.error.HTTPError or urllib.error.URLError as e:
            print(e)
        except http.client.IncompleteRead or http.client.RemoteDisconnected as e:
            if num_retries == 0:
                return
            else:
                requestImg(url, file_name, num_retries - 1)
    
    
    
    if __name__ == "__main__":
    	
    	# 基础瓦片URI
        base_url = 'http://t0.tianditu.gov.cn/img_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=img&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TILEMATRIX=11&TILEROW=%s&TILECOL=%s&tk=xxxxx'
    	
    	# 给定瓦片行列范围，构造完整瓦片URI并下载
        for r in range(893, 896):
            for c in range(1668, 1671):
                print(r, c)
                url = base_url%(r, c)
                print(url)
                file_name = "./images/" + str(r) + "-" + str(c) + '.jpg'
                requestImg(url, file_name)
                time.sleep(random.random())
    
        print("------ok------")

## 2.2 拼接卫星影像瓦片数据 ##

下载好的文件按照`ROW_COL.jpg`格式统一保存在`imagesMerge`文件夹中，需按照行列号次序将瓦片图片拼接成完整的瓦片图片。基于`Python PIL库`实现此功能。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JhaWR1XzI2NjQ2MTI5_size_16_color_FFFFFF_t_70 1]  
示例代码如下：

# -*- coding: utf-8 -*-
    #!/usr/bin/env python
    import glob, re
    from PIL import Image
    
    # 按照x、y顺序对文件名进行排序
    files = glob.glob('./imagesMerge/*.jpg')
    files.sort(key=lambda x: tuple(int(i) for i in re.findall(r'\d+', x)[:2]))
    
    # 将每一行文件名保存到一个数组中
    imagefiles = { }
    for item in files:
        match = re.search(r'\d+', item)
        pre = int(match.group())
    
        if not imagefiles.get(pre):
            imagefiles[pre] = []
            
        imagefiles[pre].append(item)
    
    # 键值对转排序后的列表
    imagefiles = sorted(zip(imagefiles.keys(), imagefiles.values()))
    
    # 预先生成合并后大小的空图片
    total_width = len(imagefiles) * 256
    total_height = len(imagefiles[0][1]) * 256
    new_image = Image.new("RGB", (total_width, total_height))
    
    # 逐行拼接
    y_offset = 0
    for item in imagefiles:
        x_offset = 0
        images = list(map(Image.open, item[1])) # 映射函数，返回列表
        for subitem in images:
            new_image.paste(subitem, (x_offset, y_offset))
            x_offset += subitem.size[0]
        y_offset += images[0].size[0]
    
    new_image.save('merge.jpg', quality = 100)

## 2.3 增加坐标信息 ##

步骤`2.2`中拼接好的图片为`JPG`格式，使用`GDAL`库完成格式转换及地理信息导入。示例代码如下：

from osgeo import gdal
    
    gcp_items = [
        [113.203125, 22.431340156360612, 0, 0],         # 左上，0行0列
        [113.73046875, 22.431340156360612, 767, 0],     # 右上，0行x列
        [113.203125, 21.943045533438177, 0, 767],       # 左下，y行0列
        [113.73046875, 21.943045533438177, 767, 767]    # 右下，y行x列
    ]
    
    gcp_list = []
    for item in gcp_items:
        x, y, pixel, line = item
        z = 0
        gcp = gdal.GCP(x, y, z, pixel, line)
        gcp_list.append(gcp)
    
    options = gdal.TranslateOptions(format='GTiff', outputSRS='EPSG:4326',GCPs=gcp_list)
    gdal.Translate('OK2.tiff', 'merge.jpg', options=options)

# 3 代码封装 #

第二节虽然实现了下载瓦片、拼接瓦片、格式转换等基本功能，但所有代码均为硬编码，缺乏交互性及灵活性。比如：

*  下载瓦片时需要人工判断感兴趣区域行列号范围
 *  拼接图片四周坐标需要手动获取

因此，需要完成感兴趣区域经纬度范围获取，经纬度及行列号转换等功能。

## 3.1 感兴趣区域获取 ##

基于`天地图Javascript API`实现矩形绘制及坐标范围获取，示例代码如下：

<!DOCTYPE html>
    <html>
    
    <head>
        <meta http-equiv="content-type" content="text/html; charset=utf-8" />
        <title>天地图影像下载</title>
        <link rel="stylesheet" type="text/css" href="/static/css/tianditu.css">
        <script type="text/javascript" src="http://api.tianditu.gov.cn/api?v=4.0&tk=abd833b69fac00096820aca1c93804ca"> </script>
        <script src="/static/js/tianditu.js"></script>
    </head>
    
    <body onLoad="onLoad()">
        <div id="mapDiv"></div>
        <div class="container">
            <div>
                <input class='bound_selector' type="button" value="框选影像范围" onClick='openRectangleTool()' />
            </div>
            <form id="coord_form">
                <input id="tl_lng" class='coord_input' type="text" name="tl_lng" placeholder="左上经度" required="required">
                <input id="tl_lat" class='coord_input' type="text" name="tl_lat" placeholder="左上纬度" required="required">
                <input id="br_lng" class='coord_input' type="text" name="br_lng" placeholder="右下经度" required="required">
                <input id="br_lat" class='coord_input' type="text" name="br_lat" placeholder="右下纬度" required="required">
                <input id="zoom" class='coord_input' type="text" name="zoom" placeholder="缩放级别" required="required" value="12">
                <input id="submit" name="submit" type="button" value="生成卫片" onClick='submitCoordForm()'>
            </form>
            <form id="result_form">
                <input id="download_file_name" type="text" name="file_name" placeholder="缩放级别" required="required" value="">
                <input class='image_downloader' type="button" value="下载卫片" onClick='downloadImage()' />
            </form>
        </div>
    </body>
    
    </html>

var map, handler;
    var zoom = 12;
    var lay;
    var onlyMapLay;
    
    
    /** * 加载地图 */
    function onLoad() { 
        var imageURL = "http://t0.tianditu.gov.cn/img_w/wmts?" +
            "SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=img&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles" +
            "&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}&tk=xxxxx";
    
        //创建自定义图层对象
        lay = new T.TileLayer(imageURL, { 
            minZoom: 1,
            maxZoom: 18
        });
        var config = { 
            layers: [lay]
        };
    
        //初始化地图对象
        map = new T.Map("mapDiv", config);
        map.centerAndZoom(new T.LngLat(113.40769, 22.20945), zoom);
        map.enableScrollWheelZoom();
    
        //创建标注工具对象
        polygonTool = new T.PolygonTool(map, config);
    }
    
    
    /** * 绘制矩形 */
    function openRectangleTool() { 
        map.clearOverLays();
        if (handler) handler.close();
        handler = new T.RectangleTool(map, { 
            follow: true
        });
        handler.open();
        handler.addEventListener("draw", onDrawRectangle);
    }
    
    
    /** * 矩形绘制完成，传递坐标 */
    function onDrawRectangle(e) { 
        var tr_lat = e.currentBounds.Lq.lat;
        var tr_lng = e.currentBounds.Lq.lng;
        var bl_lat = e.currentBounds.kq.lat;
        var bl_lng = e.currentBounds.kq.lng;
        document.getElementById("tl_lng").value = bl_lng;
        document.getElementById("tl_lat").value = tr_lat;
        document.getElementById("br_lng").value = tr_lng;
        document.getElementById("br_lat").value = bl_lat;
    }
    
    
    function submitCoordForm() { 
        
        var elements = document.getElementsByClassName("coord_input");
        var postData = (function (obj) { 
            var str = "";
            for (var i = 0; i < obj.length; i++) { 
                str += obj[i].name + "=" + obj[i].value + "&"
            }
            return str;
        })(elements);
    
        var xhr = new XMLHttpRequest();
        xhr.open("post", "/image");
        xhr.setRequestHeader("Content-type", "application/x-www-form-urlencoded");
    
        xhr.onreadystatechange = function () { 
            if (xhr.readyState == 4 && xhr.status == 200) { 
                document.getElementById("download_file_name").value = xhr.responseText;
                alert("卫片已生成！")
            }
        }
        xhr.send(postData);
    }
    
    
    function downloadImage(){ 
        var filename = document.getElementById("download_file_name").value;
        var xhr = new XMLHttpRequest();
        url = "/download/" + filename;
        xhr.open("get", url, true);
        xhr.responseType = "blob";
        xhr.setRequestHeader("Content-Type","application/x-www-form-urlencoded");
    
        xhr.onreadystatechange = function () { 
            if (xhr.readyState == 4 && xhr.status == 200) { 
                var content = xhr.response;
                var file_url = window.URL.createObjectURL(new Blob([content]), { type: ' image/tiff'})
                var link = document.createElement("a");
                link.style.display = "none";
                link.href = file_url;
                link.download = xhr.getResponseHeader('filename');
                document.body.appendChild(link);
                link.click();
                document.body.removeChild(link);
            }
        }
        xhr.send();
    }

body, html { 
        position: relative;
        width: 100%;
        height: 100%;
        margin: 0;
        font-family: "Microsoft YaHei"
    }
    
    #mapDiv { 
        width: 100%;
        height: 100%
    }
    
    .container { 
        position: absolute;
        right: 0;
        top: 0;
        z-index: 9999;
        width: 250px;
        height: 100%;
        background: rgba(255, 255, 255, 1);
    }
    
    .bound_selector, .coord_input, #submit, .image_downloader { 
        display: inline-block;
        box-sizing: border-box;
        width: calc(100% - 40px);
        height: 35px;
        margin: 5px 20px;
    }
    
    .bound_selector, #submit, .image_downloader { 
        border: none;
        background: #000;
        color: #fff;
        cursor: pointer;
    }
    
    .bound_selector:hover, #submit:hover, .image_downloader:hover { 
        background: #555;
    }
    
    .coord_input { 
        padding: 5px;
        border: 2px solid #aaa;
    }
    
    #download_file_name { 
        display: none;
    }

## 3.2 经纬度行列号转换 ##

# -*- coding: utf-8 -*-
    #!/usr/bin/env python
    import math
    
    def deg2num(lon_deg, lat_deg, zoom):
      lat_rad = math.radians(lat_deg)
      n = 2.0 ** zoom
      ytile = int((lon_deg + 180.0) / 360.0 * n)
      xtile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
      print(xtile, ytile)
    
    def num2deg(xtile, ytile, zoom):
      n = 2.0 ** zoom
      lon_deg = ytile / n * 360.0 - 180.0
      lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * xtile / n)))
      lat_deg = math.degrees(lat_rad)
      print(lon_deg, lat_deg)
    
    
    if __name__ == "__main__":
        zoom = 11
        x = 896
        y = 1671
    
        num2deg(x, y, zoom)
        deg2num(113.73046875,21.943045533438177,11)

## 3.3 Flask App ##

将上述功能模块封装成一个Flask App，[项目地址][Link 1]。

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JhaWR1XzI2NjQ2MTI5_size_16_color_FFFFFF_t_70]: /images/20221022/c202c9a98a9f47568a07a87367eb3487.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JhaWR1XzI2NjQ2MTI5_size_16_color_FFFFFF_t_70 1]: /images/20221022/dc33cee9a3174c59a5a88284d8314c92.png
[Link 1]: https://github.com/moyemoji/TiandituApp