一、原理
1.瓦片地图原理：瓦片地图原理- 简书 ()
二、过程
爬取数据
1.找到矢量瓦片服务地址，以及瓦片的请求规则，构造请求url
2.计算瓦片范围，通过查看服务参数信息，发现服务为arcgisserver服务，
服务信息地址：Services Directory - Map_4(VectorTileServer) ()
可以看到范围，瓦片层级对应的容差和比例，瓦片原点信息，利用这些信息可以用来计算瓦片的位置（我这里使用代码计算，二次验证瓦片信息的准确性）：
3.下载瓦片，下载的瓦片是pbf格式，因为坐标都在一个瓦片范围内，并且图层都冗杂在一起，需要解析。
4.解析pbf，提取Preliminary Floodplain图层矢量信息，我们知道了每个瓦片的大小，现在相对行列号最小的瓦片进行偏移坐标，使图层按正确的方式排列
重新排列后的图层
5.使用Qgis把图层先进行合并—修正几何图形(合并后的图层可能会有一些误差，需要修复一下)—融合操作(融合字段选中mvt_id)—导出shp，导出shp后需要把prj文件删除掉，不然在arcgis中无法查看
合并输出
修正几何图层
融合图层
导出图层
删除prj坐标系文件
6.空间校正，融合导出后的图层坐标系是错误的，需要进行空间校正，这里我们使用arcgis的空间校正工具进行校正，我们先在正确的图层上取出四个点来进行校正。
1）取点，首先需要建立arcgis矢量瓦片图层的连接，右键Vector Tiles—新建通用连接—输入名称Map_4,输入瓦片服务的链接/tiles/V3rkP5g6N5bDtF74/arcgis/rest/s，前面我们也提到过，将最大层级设置到21级，放大图层取四个特征点，最好处于四个角，方便空间校正。
2）在arcgis中打开图层，我们把复制的坐标创建一个点shp，方便校正，关于空间校正方法可以看下这里ArcGIS 空间校正-百度经验 ()
校正后图层
最后，虽然图层已经得到，但是仍有缺失得地方，这时我们需要重复操作以上部分步骤，获取其他级别得图层来进行补缺。
三、代码运行
1.安装gdal和requests库
2.更改代码的文件放置位置，78行改成你自己的文件夹位置，不需要存在，代码检查到不存在会创建的。
先运行download_tiles函数，运行后会在你输入的目录下生成数据目录，里面放置了pbf图层瓦片。
再运行calc_tiles_location函数，会在目录下生成解析并重新拼接后的geojson数据，该数据需要qgis打开，arcgis无法打开。
import mathimport requestsimport osfrom osgeo import ogrfrom osgeo.gdalconst import GA_ReadOnlyimport jsonimport 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 = '/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 Falsetry:# 创建文件夹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]) / resolutiony_pos = (origin[1] - coord_y) / resolutionx_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 * resolutionxmin = x * tile_size * resolution + origin[0]ymin = origin[1] - (y + 1) * tile_size * resolutionreturn xmin, ymin, xmax, ymaxdef download_tiles():# 获取瓦片的层级，我这里使用了4个层级的数据，因为面关系复杂，所以需要多个图层层级比对14，11，10，8z = 10# 爬取后发现全都在一个地图，需要重新计算瓦片位置，每个瓦片的大小是4096width = 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 - 1start_y_tile = start_y_tile + 1end_x_tile = end_x_tile + 1end_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，8z = 10# 爬取后发现全都在一个地图，需要重新计算瓦片位置，每个瓦片的大小是4096width = 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 - 1start_y_tile = start_y_tile + 1end_x_tile = end_x_tile + 1end_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) * widthy_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 sourcelayer_name = layer.GetName() # get layer namefeatures = []result_geojson_str = {"type": "FeatureCollection","features": []}if layer_name == 'Preliminary Floodplain':for feature in layer: # get feature in layerm_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'] = featureswith open(result_geojson, encoding='utf-8', mode='w') as f:json.dump(result_geojson_str, f)# 重新计算瓦片位置if __name__ == '__main__':# 先运行这个函数download_tiles()# calc_tiles_location()
