干GIS这行久了你会发现“下载瓦片再合并”是个永远躲不开的活。区县级瓦片下载与合并TIF、PNG这种需求规划、应急、测绘、政务项目里能碰到几百次领导要一张离线底图项目要一套能叠到自己数据上的影像背景或者甲方工程区压根没网。最近我做的就是这件事——把一个区县边界范围内的在线瓦片整套拉下来合并成一张带坐标的TIF和一张干净的PNG。说实话原理不复杂但瓦片清单怎么算、并发拉多少会被限流、八千张图拼进一张图内存怎么扛、PNG怎么带上坐标变成TIF每一步都是细节活。这篇把我完整跑通的方案和踩过的坑都写出来正在做离线底图或批量栅格合并的朋友可以直接抄作业。1. 需求拆解与方案选型1.1 先想清楚你究竟要下载什么、合并成什么区县级瓦片下载这事儿本质上是在回答三个问题范围多大、要几级、输出什么格式。范围层面区县边界通常是一个面要素可能是SHP、GeoJSON、KML或者数据库里的一个多边形。你要做的不是下载这个面本身而是下载能完全覆盖这个面的所有地图瓦片。瓦片是正方形的所以哪怕边界是个歪歪扭扭的区县形状实际下载范围永远是边界的外接矩形Bounding Box。区别只在于后续你是保留矩形底图还是再按行政区边界做一次裁剪。层级层面区县级一般要14到18级。14级适合看全区轮廓16级能看清道路和村落18级基本就是到单体建筑了。我这次的项目主要做16级兼顾15级总共三个层级瓦片数量已经上万。不同层级瓦片数量按4倍递增选层级前先估一下量不然下载时间和磁盘空间都会失控。输出格式层面PNG和TIF各有用处。PNG适合做汇报图、底图叠加、网页端展示文件体积小、打开快TIF更准确说是GeoTIFF自带地理坐标信息适合放进ArcMap、QGIS这类专业软件里和其他数据精确套合。很多人以为两张图内容一样随便输出一个就行实际上坐标信息和位深差异会让TIF在很多场景里不可替代。1.2 三条路径怎么选插件、图形工具、脚本实现“下载合并”有三条常见路线我分别试过各有各的坑。第一条是用QGIS这类开源桌面软件装个瓦片插件比如QuickMapServices或同类工具加载底图后按画布范围导出。这条路最快适合一次性小范围出图但问题也明显导出范围只能靠手动画框没法精确绑定到某个面要素层级多了以后QGIS会把所有瓦片重投影再合成一次速度慢得让人想关机而且QGIS内部处理透明通道的方式会自动填充背景出来的PNG经常有奇怪的底色残边。第二条是用现成的商业化地图下载工具。这类工具通常号称“输入范围自动下载、自动拼接”确实省事但我在实际项目里被坑过两次一次是坐标体系不透明工具按自己的墨卡托变体下载合并出来和本地CGCS2000坐标系数据差了十几米另一次是下载策略写死目标服务器稍微限流它就疯狂重试最后账号被临时封禁。工具能用但你把控制权交出去了出问题只能干瞪眼。第三条就是我最终采用的脚本化方案用Python控制下载、拼接、配准的每一个环节。初始写脚本要多花几个小时但换来的是范围精确可控、并发速度可调、格式输出自由、同一个区县换边界换层级直接复用。对经常处理同类需求的人来说脚本方案的综合成本反而是最低的。1.3 为什么最终选择了脚本化全链路选择脚本化不是因为“写代码更高级”而是因为瓦片处理链条里的每一个环节都需要精确控制。下载环节需要精确的并发数和重试逻辑既要跑得快又不能被服务器判定为恶意请求合并环节需要自主控制画布大小和透明通道处理策略避免QGIS那种自动填背景配准环节需要自己算出准确的像元大小和坐标原点这直接决定TIF叠加是否偏移。这些细节图形工具全都封装死了你只能用默认逻辑出了问题都不知道去哪改。脚本方案还有个意外的好处可审计。下载了哪些瓦片、多少张成功、多少张重试都有日志出了问题能定位到具体某一个z/x/y编号。这种确定性在项目交付阶段非常宝贵因为甲方问“你这张图覆盖全了吗”的时候你能拿出瓦片清单给他看而不是说“应该是全了吧”。2. 关键原理瓦片编号与边界范围计算2.1 瓦片金字塔与z/x/y编号规则地图瓦片采用金字塔结构每一级zoom记为z把全球地图切成4的z次方个256×256的小方块。x是从左到右的列号y是从上到下的行号。这套规则所有主流在线地图服务都差不多但坐标原点和投影方式有差异后面细说。既然要按区县边界下载第一步就是把经纬度范围换算成瓦片编号范围。标准Web墨卡托EPSG:3857的换算公式是固定的核心就两句import math def lonlat_to_tile(lon, lat, z): n 2 ** z x int((lon 180.0) / 360.0 * n) # y方向用的是墨卡托投影后的纬度不是直接线性换算 lat_rad math.radians(lat) y int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n) return x, y注意y坐标不能直接用“纬度占比乘n”那样在高纬度地区会严重变形。必须先把纬度做Web墨卡托正算再换算成行号。很多人在这一步写错导致下载区域南北方向整体偏移合并出来的图和高精度矢量数据完全对不上。反向换算也不难用来把瓦片编号变回经纬度范围合并时写坐标信息要用def tile_to_lonlat(x, y, z): n 2 ** z lon_left x / n * 360.0 - 180.0 lat_rad math.atan(math.sinh(math.pi * (1.0 - 2.0 * y / n))) lat_top math.degrees(lat_rad) lon_right (x 1) / n * 360.0 - 180.0 lat_rad math.atan(math.sinh(math.pi * (1.0 - 2.0 * (y 1) / n))) lat_bottom math.degrees(lat_rad) return lon_left, lat_top, lon_right, lat_bottom2.2 从区县边界矢量算出瓦片下载清单有了换算函数剩下的就是纯计算。先读取边界矢量文件拿到外接矩形经纬度范围再对每个目标层级分别算出瓦片行列号范围。def calc_tile_range(bbox, z): lon_min, lat_min, lon_max, lat_max bbox x_min, y_max lonlat_to_tile(lon_min, lat_max, z) x_max, y_min lonlat_to_tile(lon_max, lat_min, z) return (x_min, y_min, x_max, y_max)这里有个细节外接矩形的左下角对应的是y最大的瓦片右上角对应y最小的瓦片所以算出来要再规整一次。实际下载时遍历tiles [] for z in zoom_levels: x_min, y_min, x_max, y_max calc_tile_range(bbox, z) for x in range(x_min, x_max 1): for y in range(y_min, y_max 1): tiles.append((z, x, y))区县级范围在16级通常几十到一百多张宽三个层级加起来小两万张瓦片下载时间取决于并发和服务器响应速度一般十几分钟到一两个小时。2.3 边界裁剪矩形底图还是严格行政边界下载范围是外接矩形合并出来自然也带矩形外框。如果你的交付要求“底图必须贴合区县边界、边界外不能有内容”就得在合并之后再叠一次裁剪。ArcMap中实现这个有两种常见操作栅格裁剪Clip和按掩膜提取Extract by Mask。很多人问这俩到底啥区别我在这里顺便说清楚Clip是数据管理工具箱里的工具按面要素范围裁剪栅格如果勾选了“使用输入要素裁剪几何”输出范围会贴合面形状但输出栅格仍然按外包矩形存储Extract by Mask是空间分析工具箱的工具它对整个栅格逐像元判断面外的像元直接写成NoData输出栅格理论上仍是矩形但面外都是无效值叠加显示时只显示面内部分。实际做底图交付我一般先用脚本合并出矩形底图再交给内业人员在ArcMap里用Extract by Mask做精确裁剪。这样边界处是真正的NoData后续做统计分析不会把边界外像元算进去。注意两种工具的像元对齐策略不同裁剪前先确认参考栅格避免输出分辨率被悄悄改掉。2.4 不同底图源的坐标体系差异这一步是坑最多的地方。OSM、标准Web墨卡托底图用的是WGS84地理坐标对应的EPSG:3857投影天地图虽然也用Web墨卡托但坐标基准是CGCS2000百度地图的瓦片在Web墨卡托基础上又做了偏移加密处理而且百度瓦片编号的z值定义和谷歌系不一样缩放级别差一级。换句话说同一个经纬度点在不同服务里算出的瓦片编号是不同的一套数字。所以下载前必须确认你用的瓦片源到底遵循哪套规则瓦片源坐标系瓦片z定义合并后TIF投影OSM/标准XYZWGS84(Web墨卡托)与全球金字塔一致EPSG:3857天地图CGCS2000(Web墨卡托)与全球金字塔基本一致CGCS2000 / EPSG:3857百度地图BD09坐标(加密偏移)起始级别有偏移需按百度规则写高德/腾讯GCJ02坐标(加密偏移)与全球金字塔基本一致GCJ02对应投影我的建议是优先选择标准XYZ服务作为底图源因为后续投影、坐标转换最省心。如果项目必须在特定平台上取图那下载脚本里的换算公式、输出TIF的投影定义都要按那个服务的规则单独写一套不能拿通用公式硬套。3. 实操环节从并发下载到合并输出3.1 环境准备与依赖安装我用的Python 3.10 四个库requests负责下载Pillow处理图像拼接GDAL写GeoTIFFtqdm显示进度。安装命令很简单pip install requests pillow tqdm gdalGDAL在Windows上偶尔装不顺利装不上就换conda装conda install gdal或者用rasterio替代功能上够用。另外推荐装一个shapely后面如果要按面边界过滤瓦片中心点用它做空间判断很方便。3.2 并发下载速度与限流的平衡下载瓦片最忌讳两件事一点点串行下下到猴年马月或者并发开满把服务器打爆然后被限流。我实测下来的经验是目标服务器一般容忍4到8个并发超过10个很容易触发限流。脚本里固定线程池大小5每条线程独立拉取失败自动重试。from concurrent.futures import ThreadPoolExecutor, as_completed import requests, time HEADERS {User-Agent: Mozilla/5.0 (GIS project; tile downloader)} def download_tile(url, local_path, retries3): for attempt in range(retries): try: resp requests.get(url, headersHEADERS, timeout10) if resp.status_code 200 and len(resp.content) 100: with open(local_path, wb) as f: f.write(resp.content) return True elif resp.status_code 429 or resp.status_code 503: time.sleep(2 * (attempt 1)) except requests.RequestException: time.sleep(2) return False def download_all(tile_list, url_template, out_dir, workers5): tasks [] with ThreadPoolExecutor(max_workersworkers) as pool: for z, x, y in tile_list: local_path f{out_dir}/{z}/{x}/{y}.png url url_template.format(zz, xx, yy) tasks.append(pool.submit(download_tile, url, local_path)) for future in as_completed(tasks): if not future.result(): print(下载失败瓦片:, future)三个细节分享一是必须带User-Agent很多服务器会拦截空请求头二是判断成功的条件不只要看状态码还要看响应体大小有些服务对不存在的瓦片会返回一张200的空图或占位图体积只有几百字节这种要视为失败三是有条件的话在两次请求之间加一个极短的随机延迟比如0.1到0.3秒能显著降低被限流概率。3.3 合并PNG画布拼接与透明通道处理所有瓦片下齐后合并就是纯粹的图像操作。先按瓦片数量创建画布然后按行列号把每张瓦片粘贴到对应位置from PIL import Image def merge_tiles(tiles_dir, z, x_min, y_min, x_max, y_max, out_png): width (x_max - x_min 1) * 256 height (y_max - y_min 1) * 256 canvas Image.new(RGBA, (width, height), (0, 0, 0, 0)) for x in range(x_min, x_max 1): for y in range(y_min, y_max 1): tile_path f{tiles_dir}/{z}/{x}/{y}.png try: tile Image.open(tile_path).convert(RGBA) except FileNotFoundError: continue canvas.paste(tile, ((x - x_min) * 256, (y - y_min) * 256)) canvas.save(out_png, PNG)注意几个坑第一画布模式务必用RGBA而不是RGB因为很多底图瓦片本身带透明通道RGB模式会把透明区域变成黑色或者白色残边第二paste时如果瓦片本身是RGBA默认会直接把原始RGBA粘上去不会做alpha混合对纯色底图没问题但如果瓦片边缘有抗锯齿渐变连续拼接处可能出现细线这时候改用canvas.alpha_composite(tile, position)会过渡得更自然代价是性能下降第三遇到缺图不要直接崩溃先跳过最后统一补下。3.4 从PNG到TIF坐标配准才是关键PNG本身不带地理坐标要让它在ArcMap里精确叠加就得写一个世界文件PNG对应的扩展名是.pgw或者直接输出GeoTIFF。GeoTIFF更省事我建议直接写TIF。核心是算对像元大小和坐标原点。Web墨卡托全球东西跨度为40075016.68557849米某一层级下单个像元对应地面大小为pixel_size 40075016.68557849 / (256 * (2 ** z))合并图左上角的坐标就是该层级第x_min列瓦片左边缘和第y_min行瓦片上边缘对应的投影坐标import math from osgeo import gdal, osr def merge_to_geotiff(png_path, z, x_min, y_min, x_max, y_max, out_tif): src Image.open(png_path) pixel_size 40075016.68557849 / (256 * (2 ** z)) origin_x -20037508.342789244 x_min * 256 * pixel_size origin_y 20037508.342789244 - y_min * 256 * pixel_size driver gdal.GetDriverByName(GTiff) ds driver.Create(out_tif, src.width, src.height, 4, gdal.GDT_Byte) ds.SetGeoTransform([origin_x, pixel_size, 0, origin_y, 0, -pixel_size]) srs osr.SpatialReference() srs.ImportFromEPSG(3857) ds.SetProjection(srs.ExportToWkt()) for band in range(4): data src.getchannel(band).tobytes() ds.GetRasterBand(band 1).WriteRaster(0, 0, src.width, src.height, data) ds.FlushCache()这里最容易被忽略的是“原点”的定义。GDAL的GeoTransform里第0、3个参数分别表示栅格左上角的投影X、Y坐标不是栅格中心坐标。默认情况下瓦片栅格的原点就是瓦片左上角你算出的origin_x、origin_y直接对应合并图左上角第一张瓦片的左上角完全没差。但如果你拿某个地图服务的元数据直接抄人家给的可能就是中心点坐标那叠加时就会偏半个像元图上看着不明显量距离误差却实实在在。3.5 大范围合并内存爆炸的解决方案上面的直接合并方式有个天花板图片太大时内存扛不住。一张16级、横竖100张瓦片的图就是25600×25600像素RGBA模式占内存接近2.6GB普通办公电脑直接卡死。解决思路是把“一次性拼大图”改成“分块拼、再拼接”。最简单的一个做法是把大范围按四象限切成四块分别合并成四张子图再统一贴到最终画布上。更工程化的做法是给每一张瓦片单独生成带坐标的小GeoTIFF然后用gdalbuildvrt把几百上千个小文件构建成虚拟栅格最后用gdal_translate落盘成一张大的Tiled TIF。这个方案内存占用极其稳定因为GDAL内部做了分块读取不会把整张大图一次性load进内存。# 先把每个瓦片配准成小GeoTIFF脚本循环处理 # 然后用虚拟栅格合并 gdalbuildvrt merged.vrt tile_mosaic_dir/*.tif gdal_translate -co TILEDYES -co COMPRESSDEFLATE merged.vrt final.tif用这种方式合并出来的TIF是分块存储的ArcMap读取效率也非常高。如果你要的层级更高、范围更大我建议从一开始就放弃PIL大图方案直接用GDAL这套。4. 常见问题与排查技巧实录4.1 下载中途大量失败或被限流最常见的表现是下载到一半突然几百张连续失败或者请求响应时间明显变长。这基本就是并发太大或者请求频率太高导致的。排查方法看失败瓦片是否是连续的x、y区间如果是基本确定为限流如果失败区域随机分布多半是网络波动或单张瓦片服务器偶发错误。处理手法就两条限制并发到3到4重试间隔翻倍。我用的是指数退避第一次等2秒、第二次4秒、第三次8秒三次还失败就把瓦片号写进failed.txt所有下载任务结束后单独跑一轮重试。重试那一轮并发降到2成功率能回到99%以上。千万别在限流期间硬刚等几十分钟再补效果最好。4.2 合并后出现白色缝隙或黑色边缘这个问题我排查过很多次原因有三种。第一种是瓦片本身带透明边缘无缝拼接时透明区域叠不出颜色远看像白色细线。处理方式是把画布背景先填充成浅灰色或瓦片源底色再用alpha_composite叠加基本能消除。第二种是某些瓦片下载不完整PNG只有上半部分有内容下半部分是黑的或空的。这种多半是下载时请求被中断但状态码返回了200被我的“响应体大于100字节”的判断漏过去了。解决方法是下载完成后做一个完整性校验每张PNG解码后的尺寸必须是256×256解析异常的直接重下。第三种是瓦片源本身就缺图尤其是偏远的区县边缘个别层级存在空白瓦片。这种没法靠重试解决只能接受或者降低一个层级取图来补。4.3 画布尺寸超限导致PIL报错Pillow在Windows上对超大图像的尺寸有限制超过一定像素尺寸会直接抛异常。我在处理18级全县范围时就遇到过报错信息大意是“image size exceeds limit”。解决方式有两个一个是用Image.MAX_IMAGE_PIXELS None关掉限制但这不是根治内存照样爆另一个就是前面说的分块合并或GDAL VRT方案。我个人强烈建议直接上GDAL路线因为内存可控是第一位的PIL大画布方案只适合几百张瓦片的小范围。4.4 TIF叠加后整体偏移半张瓦片这是坐标配准最容易出的问题而且特别隐蔽。表现是TIF叠加到ArcMap里和矢量边界整体错位放大看大概偏移半个瓦片到一整个瓦片的距离而且越放大偏差越大量出来是固定米数而不是固定角度。原因基本就两个一是GeoTransform原点坐标写错用了瓦片中心而不是左上角二是投影定义错了比如瓦片源是GCJ02坐标你却写了EPSG:3857两个坐标系本身就有固定的偏移。排查步骤很直接先用gdalinfo看TIF的坐标范围是否和边界矢量范围吻合再用一个已知地物点核对图上位置最后确认输出投影和底图源是否同坐标系。我自己踩过最大的坑是拿天地图瓦片合并后直接标成EPSG:3857忘记了天地图基准是CGCS2000和WGS84在区县级范围虽然差异不大但叠加高精度测量数据时还是能看出几米的偏差。4.5 常见问题速查表现象可能原因解决方法下载中途大量失败并发过高被限流降并发、指数退避、错峰重试响应码200但文件损坏请求中断、占位图校验解码尺寸异常重下合并后有白色细线透明通道处理不当画布填充底色 alpha_composite合并后有大片黑色瓦片源缺图或404换层级、找替代源补图PIl直接报size exceeded画布超限改GDAL VRT分块方案TIF叠加整体偏移原点坐标或投影写错核对GeoTransform和投影定义边界外显示空白未做掩膜提取转GeoTIFF后用Extract by Mask输出TIF被拉伸变形像元大小算错重新按公式计算pixel_size5. 经验之外的几个提醒最后想分享两个实操中总结出的经验不算技术但很影响项目体验。第一下载任何一种在线地图瓦片前先想清楚使用场景和版权边界。内部项目做底图分析一般问题不大但如果要对外发布或者做商业产品不同地图服务的授权条款差别很大。我的习惯是内部项目优先选开源或政务类底图对外项目提前和法务确认授权别等图都合并完了才想起来版权的事。第二脚本要写得“可续跑”。下载两万张瓦片过程中断网、断电、电脑重启都可能发生。所以我的脚本每次运行前会扫描本地已有的瓦片文件跳过已存在的只补缺失的。这个逻辑加上前面提到的failed.txt补充重试机制能让整个下载过程随时中断、随时继续不浪费任何一次网络请求。这个习惯救过我很多次毕竟两万张瓦片重新拉一遍谁顶得住。这个流程本质上就是一个可复用的离线底图生产管线输入一个面边界和层级列表输出一个带坐标的TIF和一个干净的PNG。下次接到类似需求改个边界、调个层级、换个瓦片源URL脚本就能直接再跑一遍。工具化的好处就在这里一次投入后面全是复用。