|
一、原理
1.瓦片地图原理：瓦片地图原理- 简书 (jianshu.com)
二、过程
爬取数据
1.找到矢量瓦片服务地址，以及瓦片的请求规则，构造请求url
2.计算瓦片范围，通过查看服务参数信息，发现服务为arcgisserver服务，
服务信息地址：Services Directory - Map_4(VectorTileServer) (arcgis.com) 
可以看到范围，瓦片层级对应的容差和比例，瓦片原点信息，利用这些信息可以用来计算瓦片的位置（我这里使用代码计算，二次验证瓦片信息的准确性）：
3.下载瓦片，下载的瓦片是pbf格式，因为坐标都在一个瓦片范围内，并且图层都冗杂在一起，需要解析。
4.解析pbf，提取Preliminary Floodplain图层矢量信息，我们知道了每个瓦片的大小，现在相对行列号最小的瓦片进行偏移坐标，使图层按正确的方式排列
重新排列后的图层
5.使用Qgis把图层先进行合并—修正几何图形(合并后的图层可能会有一些误差，需要修复一下)—融合操作(融合字段选中mvt_id)—导出shp，导出shp后需要把prj文件删除掉，不然在arcgis中无法查看
合并输出 
修正几何图层
 
融合图层
导出图层
删除prj坐标系文件
6.空间校正，融合导出后的图层坐标系是错误的，需要进行空间校正，这里我们使用arcgis的空间校正工具进行校正，我们先在正确的图层上取出四个点来进行校正。
1）取点，首先需要建立arcgis矢量瓦片图层的连接，右键Vector Tiles—新建通用连接—输入名称Map_4,输入瓦片服务的链接https://tiles.arcgis.com/tiles/V3rkP5g6N5bDtF74/arcgis/rest/s，前面我们也提到过，将最大层级设置到21级，放大图层取四个特征点，最好处于四个角，方便空间校正。
2）在arcgis中打开图层，我们把复制的坐标创建一个点shp，方便校正，关于空间校正方法可以看下这里ArcGIS 空间校正-百度经验 (baidu.com)
 
校正后图层
最后，虽然图层已经得到，但是仍有缺失得地方，这时我们需要重复操作以上部分步骤，获取其他级别得图层来进行补缺。
三、代码运行
1.安装gdal和requests库
2.更改代码的文件放置位置，78行改成你自己的文件夹位置，不需要存在，代码检查到不存在会创建的。
先运行download_tiles函数，运行后会在你输入的目录下生成数据目录，里面放置了pbf图层瓦片。
 
再运行calc_tiles_location函数，会在目录下生成解析并重新拼接后的geojson数据，该数据需要qgis打开，arcgis无法打开。
import math
import requests
import os
from osgeo import ogr
from osgeo.gdalconst import GA_ReadOnly
import json
import numpy as np
# 获得读取
driver = ogr.GetDriverByName("MVT")
driver_geojson = ogr.GetDriverByName('GeoJSON')
info_json = {"currentVersion": 10.81, "name": "Map_4", "capabilities": "TilesOnly", "type": "indexedVector",
"defaultStyles": "resources/styles", "tiles": ["tile/{z}/{y}/{x}.pbf"], "exportTilesAllowed": 'false',
"initialExtent": {"xmin": -1.0551344062892381E7, "ymin": 3422255.6049980605, "xmax": -1.0424577025847457E7,
"ymax": 3578891.2062397744, "spatialReference": {"wkid": 102100, "latestWkid": 3857}},
"fullExtent": {"xmin": -1.0551344062892381E7, "ymin": 3422255.6049980605, "xmax": -1.0424577025847457E7,
"ymax": 3578891.2062397744, "spatialReference": {"wkid": 102100, "latestWkid": 3857}},
"minScale": 2.958287637957775E8, "maxScale": 35.26553675,
"tileInfo": {"rows": 512, "cols": 512, "dpi": 96, "format": "pbf",
"origin": {"x": -2.0037508342787E7, "y": 2.0037508342787E7},
"spatialReference": {"wkid": 102100, "latestWkid": 3857},
"lods": [{"level": 0, "resolution": 78271.516964, "scale": 2.958287637957775E8},
{"level": 1, "resolution": 39135.75848199995, "scale": 1.479143818978885E8},
{"level": 2, "resolution": 19567.87924100005, "scale": 7.39571909489445E7},
{"level": 3, "resolution": 9783.93962049995, "scale": 3.6978595474472E7},
{"level": 4, "resolution": 4891.96981024998, "scale": 1.8489297737236E7},
{"level": 5, "resolution": 2445.98490512499, "scale": 9244648.868618},
{"level": 6, "resolution": 1222.992452562495, "scale": 4622324.434309},
{"level": 7, "resolution": 611.496226281245, "scale": 2311162.2171545},
{"level": 8, "resolution": 305.74811314069, "scale": 1155581.1085775},
{"level": 9, "resolution": 152.874056570279, "scale": 577790.5542885},
{"level": 10, "resolution": 76.4370282852055, "scale": 288895.2771445},
{"level": 11, "resolution": 38.2185141425366, "scale": 144447.638572},
{"level": 12, "resolution": 19.1092570712683, "scale": 72223.819286},
{"level": 13, "resolution": 9.55462853563415, "scale": 36111.909643},
{"level": 14, "resolution": 4.777314267817075, "scale": 18055.9548215},
{"level": 15, "resolution": 2.388657133974685, "scale": 9027.977411},
{"level": 16, "resolution": 1.19432856698734, "scale": 4513.9887055},
{"level": 17, "resolution": 0.597164283427525, "scale": 2256.9943525},
{"level": 18, "resolution": 0.2985821417799085, "scale": 1128.4971765},
{"level": 19, "resolution": 0.1492910708238085, "scale": 564.248588},
{"level": 20, "resolution": 0.07464553541190416, "scale": 282.124294},
{"level": 21, "resolution": 0.03732276770595208, "scale": 141.062147},
{"level": 22, "resolution": 0.01866138385297604, "scale": 70.5310735},
{"level": 23, "resolution": 0.00933069192648802, "scale": 35.26553675}]},
"maxzoom": 18,
"minLOD": 0,
"maxLOD": 15,
"resourceInfo": {
"styleVersion": 8,
"tileCompression": "gzip",
"cacheInfo": {
"storageInfo": {
"packetSize": 128,
"storageFormat": "compactV2"
}
}
},
"serviceItemId": "56d7d941b7af43069f1384e67fb8ee87", "maxExportTilesCount": 100000}
root = r'D:/projects/0705/data1'
# 下载瓦片数据
def download(z, x, y):
# 矢量瓦片服务地址
url = 'https://tiles.arcgis.com/tiles/V3rkP5g6N5bDtF74/' \
'arcgis/rest/services/Map_4/VectorTileServer/tile/{0}/{1}/{2}.pbf'.format(z, y, x)
# 请求数据
req = requests.get(url)
filename = str(x) + '_' + str(y) + '.pbf'
print('x:{0},y:{1},z{2}'.format(x, y, z))
if req.status_code != 200:
print('下载异常')
return False
try:
# 创建文件夹
path = root + r"{0}".format(z)
isExist = os.path.exists(path)
if not isExist:
os.makedirs(path)
with open(path + '/' + filename, 'wb+') as f:
f.write(req.content)
print('下载成功')
except Exception as e:
print(e)
# 图层的范围，用来计算瓦片的起始
extent = [-1.0551344062892381E7, 3422255.6049980605, -1.0424577025847457E7, 3578891.2062397744]
# 通过坐标和缩放层级来获取矢量瓦片文件
def get_tile_by_coord(coord_x, coord_y, zoom):
# 瓦片原点
origin = [-2.0037508342787E7, 2.0037508342787E7]
# 瓦片大小
tile_size = 512
# 第一级容差
min_res = (origin[1] - origin[0]) / tile_size
# 第zoom级容差
resolution = min_res / pow(2, zoom)
x_pos = (coord_x - origin[0]) / resolution
y_pos = (origin[1] - coord_y) / resolution
x_tile = int(math.ceil(x_pos / tile_size))
y_tile = int(math.ceil(y_pos / tile_size))
return x_tile, y_tile
# 通过瓦片的行列号和层级来还原瓦片的范围坐标
def get_coord_by_tile(x, y, zoom):
# 瓦片原点
origin = [-2.0037508342787E7, 2.0037508342787E7]
# 瓦片大小
tile_size = 512
# 第一级容差
min_res = (origin[1] - origin[0]) / tile_size
# 第zoom级容差
resolution = min_res / pow(2, zoom)
xmax = (x + 1) * tile_size * resolution + origin[0]
ymax = origin[1] - y * tile_size * resolution
xmin = x * tile_size * resolution + origin[0]
ymin = origin[1] - (y + 1) * tile_size * resolution
return xmin, ymin, xmax, ymax
def download_tiles():
# 获取瓦片的层级，我这里使用了4个层级的数据，因为面关系复杂，所以需要多个图层层级比对14，11，10，8
z = 10
# 爬取后发现全都在一个地图，需要重新计算瓦片位置，每个瓦片的大小是4096
width = 4096
# 计算瓦片起始范围
start_x_tile, start_y_tile = get_tile_by_coord(extent[0], extent[1], z)
end_x_tile, end_y_tile = get_tile_by_coord(extent[2], extent[3], z)
# 多加一个瓦片范围防止爬取缺失
start_x_tile = start_x_tile - 1
start_y_tile = start_y_tile + 1
end_x_tile = end_x_tile + 1
end_y_tile = end_y_tile - 1
# 开始爬取
for x in range(start_x_tile, end_x_tile):
for y in range(end_y_tile, start_y_tile):
# region 爬取瓦片
download(z, x, y)
# region 爬取瓦片
def calc_tiles_location():
extent = [-1.0551344062892381E7, 3422255.6049980605, -1.0424577025847457E7, 3578891.2062397744]
# 获取瓦片的层级，我这里使用了4个层级的数据，因为面关系复杂，所以需要多个图层层级比对14，11，10，8
z = 10
# 爬取后发现全都在一个地图，需要重新计算瓦片位置，每个瓦片的大小是4096
width = 4096
# 计算瓦片起始范围
start_x_tile, start_y_tile = get_tile_by_coord(extent[0], extent[1], z)
end_x_tile, end_y_tile = get_tile_by_coord(extent[2], extent[3], z)
# 多加一个范围防止爬取缺失
start_x_tile = start_x_tile - 1
start_y_tile = start_y_tile + 1
end_x_tile = end_x_tile + 1
end_y_tile = end_y_tile - 1
# 开始爬取
for x in range(start_x_tile, end_x_tile):
for y in range(end_y_tile, start_y_tile):
# 重新计算瓦片位置
x_diff = (x - start_x_tile) * width
y_diff = (y - end_y_tile) * width
# 放置瓦片的位置
pbf_file = root + '{2}/{0}_{1}.pbf'.format(x, y, z)
result_geojson = root + '_geojson{0}'.format(z) + '/{0}_{1}.geojson'.format(x, y)
if os.path.exists(root + '_geojson{0}'.format(z)) is not True:
os.makedirs(root + '_geojson{0}'.format(z))
if os.path.exists(pbf_file) and os.path.getsize(pbf_file) > 0:
# extent = get_coord_by_tile(x, y, z)
data_source = driver.Open(pbf_file, GA_ReadOnly)
layer_num = data_source.GetLayerCount()
for i in range(layer_num):
layer = data_source.GetLayer(i) # get layers in data source
layer_name = layer.GetName() # get layer name
features = []
result_geojson_str = {
"type": "FeatureCollection",
"features": [
]
}
if layer_name == 'Preliminary Floodplain':
for feature in layer: # get feature in layer
m_feature = json.loads(feature.ExportToJson())
geo_type = m_feature['geometry']['type']
coords = []
print(geo_type)
if geo_type == 'Polygon':
coords = m_feature['geometry']['coordinates'][0]
for c in range(len(coords)):
coords[c] = [coords[c][0] + x_diff, coords[c][1] - y_diff]
m_feature['geometry']['coordinates'] = [coords]
else:
coords = m_feature['geometry']['coordinates'][0][0]
for c in range(len(coords)):
coords[c] = [coords[c][0] + x_diff, coords[c][1] - y_diff]
m_feature['geometry']['coordinates'] = [[coords]]
features.append(m_feature)
result_geojson_str['features'] = features
with open(result_geojson, encoding='utf-8', mode='w') as f:
json.dump(result_geojson_str, f)
# 重新计算瓦片位置
if __name__ == '__main__':
# 先运行这个函数
download_tiles()
# calc_tiles_location()
|