大型TIFF 影像转MBTiles学习记录，ArcGIS pro 、ArcPy以及QGIS结合使用

先生余枫

已于 2024-06-25 10:09:50 修改

阅读量1.3k

点赞数 7

分类专栏： ArcPy ArcGIS MBTiles 文章标签：学习 arcgis python

于 2024-06-22 10:57:45 首次发布

本文链接：https://blog.csdn.net/qq_51165365/article/details/139878535

版权

ArcPy 同时被 3 个专栏收录

1 篇文章

订阅专栏

ArcGIS

1 篇文章

订阅专栏

MBTiles

1 篇文章

订阅专栏

一、数据：

1.多个(40多)位置偏移（统一偏移,偏移很大，南北半球跨度）的tif,坐标系WGS1984，原始数据很大，压缩以后15G，最后会有70多G，有多余的需要先保存其他。

2，配准图、行政区划shp各一幅，坐标系CGCS2000

二、要求：

行政区划下一级的各个地区的MBTiles，精度可观。

三、路线：

闭坑：

1.融合镶嵌 $\rightarrow$ 2. 坐标转化与投影 $\rightarrow$ 3.地理空间配准 $\rightarrow$ 4.裁剪tif $\rightarrow$ 5.转化MBTiles

修改：

第一步不应该把数据镶嵌，投影大数据（可能超过10G就慢，我20G，后面成60G）会很慢，应该先投影，再镶嵌，再校准。

1.转化与投影（批量） $\rightarrow$ 2.融合镶嵌 $\rightarrow$ 3.空间配准 $\rightarrow$ 4.裁剪tif $\rightarrow$ 5.转化MBTiles

四、工具与软件：

1.ArcGIS Pro （本人2.8或3.1版本都可行）:

1-4 步(对他熟悉一些,QGIS不熟悉，但也能做应该，对大量数据的处理，需求高）

2.Pycharm+ArcPy:

Pycharm新建工程文件，项目文件设置解释器选择pro的 python.exe ,路径：

 ArcGIS\Pro\bin\Python\envs\arcgispro-py3

有conda环境的本质上不影响，我自己的没问题，别人的我遇到了一些问题，那就是python环境无法激活，后面没有搞，直接脚本执行。

QGIS:

必须使用自带的切片【生成XYZ瓦片（MBtiles）】工具，Qtiles不建议使用，有bug，至少3.6版本不可以转化MBtiles，

建议：大量的数据对电脑性能要求高，请注意处理环境，如果是小数据则无所谓。

2.执行过程：

注意栅格源数据栅格信息要尽量保证一致，不能无脑默认。

1.栅格坐标转化与投影（ArcPy批量实现）

我需要由WGS1984到CGCS2000-3D-GK-Z-39,我使用了ArcPy,这里使用投影栅格批处理也是一样。我直接写死输入输出。

Arcpy代码坐标转化投影如下

import arcpy
import os
from arcpy import env
# 设置工作空间和环境设置
env.workspace = r"J:\TIF_res\oldTIF1984"
out_workspace = r"J:\TIF_res\your.gdb"
# 定义源坐标系统和目标坐标系统
source_coordinate_system = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;0.001;IsHighPrecision"
target_coordinate_system = "PROJCS['CGCS2000_3_Degree_GK_Zone_39',GEOGCS['GCS_China_Geodetic_Coordinate_System_2000',DATUM['D_China_2000',SPHEROID['CGCS2000',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Gauss_Kruger'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',117.0],PARAMETER['Scale_Factor',1.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"
# 遍历文件夹中的所有.tif文件
for raster_file in arcpy.ListRasters("*.tif"):
    print("投影。"+raster_file[:-4])
    # 使用 replace 函数去除所有的 '-' 符号
    sanitized_name = raster_file.replace("-", "")
    # 构建输出栅格的名称
    output_raster_name = "c" + sanitized_name[:-4] + "_ProjectRaster"
    output_raster_path = os.path.join(out_workspace, output_raster_name)
    
    # 执行投影栅格操作
    arcpy.management.ProjectRaster(in_raster=raster_file,
                               out_raster=output_raster_path,
                               out_coor_system=target_coordinate_system,
                               resampling_type="NEAREST",
                               in_coor_system=source_coordinate_system)
print("投影完成。")

2.镶嵌栅格数据至新栅格(ArcPy)

请注意mosaic_properties，参数，需要与你原来栅格一致

import arcpy
import os
import time
# 设置工作空间和输出路径
input_folder = r'J:\TIF_res\your.gdb'
output_raster_path = r'J:\TIF_res\MosaicRaster.tif'
# 定义镶嵌参数
mosaic_properties = {
    "PIXEL_TYPE": "8_BIT_UNSIGNED",  
    "NUMBER_OF_BANDS": 3,            
    "COLORSPACE_MATCHING": "FIRST",   
}
# 设置工作空间
arcpy.env.workspace = input_folder
print("处理中...")
arcpy.AddMessage("处理中...")
start_time = time.time()
# arcpy.ListRasters("*", "TIF")获取输入文件夹内所有.tif文件的路径列表 gdb应该直接使用arcpy.ListRasters()而不指定格式，所有栅格数据集
input_rasters = arcpy.ListRasters()
if input_rasters:
    # 将列表转换为字符串，方便在消息中显示
    raster_list_str = ", ".join(input_rasters)
    arcpy.AddMessage(f"找到的栅格数据集：{raster_list_str}")
    # 使用arcpy的MosaicToNewRaster函数
    arcpy.MosaicToNewRaster_management(
        input_rasters=input_rasters,
        output_location=os.path.dirname(output_raster_path),
        raster_dataset_name_with_extension=os.path.basename(output_raster_path),
        pixel_type=mosaic_properties["PIXEL_TYPE"],
        number_of_bands=mosaic_properties["NUMBER_OF_BANDS"],
        mosaic_method="LAST",  
        mosaic_colormap_mode=mosaic_properties["COLORSPACE_MATCHING"],
    )
else:
    arcpy.AddMessage("未找到任何栅格数据集。请检查input_gdb路径是否正确。")
end_time = time.time()
time_spent = end_time - start_time
arcpy.AddMessage(f"花费{time_spent}秒镶嵌完成")
print(f"任务花费了 {time_spent} 秒")
print("镶嵌完成")

3.空间配准

前提，保证配准图与需要纠正的栅格坐标系一致，然后添加控制点，或者导入你的控制点文件。

对于像我那种偏差大的，分几步来：

1> 加载号配准图与偏移图，建议配准图在最底层，利用好书签工具，

2> 取消地理配准自动应用，顺时针方向四边形形状刺点，ABCD

偏移图层源（起点） $\rightarrow$ （配准点）终点。四个方向（不一定矩形），先一阶多项式校正，应用，如果精确度还行，直接保存，偏移图层就会靠近配准图，刺点很方便了。然后慢慢增加控制点，在换另外的校正方法，自己看情况调整。

4.shp掩膜提取栅格裁剪tif（Arcpy）

这里应该是按shp属性提取按掩膜栅格(我提前把shp手动按属性分割了)，代码如下：

import arcpy
import os
# 设置工作空间和环境设置 shp文件位置
arcpy.env.workspace = "I:\TIF_res\xz_xzq"
arcpy.env.overwriteOutput = True
# 指定栅格数据集路径
raster_path = "I:\TIF_res\\MosaicRaster.tif"
# 设置输出目录
output_folder = "I:\TIF_res\\res_tif_xzq"
print("开始执行...")
arcpy.AddMessage("开始执行...")
# 遍历文件夹中的所有Shapefile
for shp_file in arcpy.ListFeatureClasses():
    # 获取乡镇名称，假设文件名为“乡镇名称.shp”，则提取“乡镇名称”
    town_name = os.path.basename(shp_file).split('.')[0]
    
    # 构建输出栅格的完整路径
    out_raster_path = os.path.join(output_folder, f"Extract_{town_name}.tif")
    
    # 确保输出目录存在
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)
    
    arcpy.AddMessage("正在执行处理" + town_name + "...")
    
    # 在处理每个shapefile之前，设置环境的处理范围为该shapefile的范围
    arcpy.env.extent = arcpy.Describe(shp_file).extent
    arcpy.env.mask = shp_file
    
    try:
        arcpy.sa.ExtractByMask(raster_path, shp_file).save(out_raster_path)
        arcpy.AddMessage(f"{town_name}按掩膜提取完成。")
    except Exception as e:
        arcpy.AddMessage(f"处理{town_name}时出错：{e}")
        
print("所有指定乡镇的按掩膜提取完成！")
arcpy.AddMessage("所有指定乡镇的按掩膜提取完成！")

tif 黑色背景掩膜不成功处理,也就是需要把栅格外的黑色，设置为noDATA,我使用手动处理会把栅格内部的黑色值夜处理，不符合需求，应该需要栅格计算器，但是我懒得用,下面是tif黑色背景掩膜代码。

import arcpy
import os
# 设置工作空间和环境设置 shp与tif一一对应关系可用
arcpy.env.workspace = "I:\\ganxian_xz_xzq"
arcpy.env.overwriteOutput = True
# 指定栅格数据集所在文件夹
raster_folder = "I:\\TIF_res\\res_tif_xzq"
raster_files = [f for f in os.listdir(raster_folder) if f.endswith('.tif')]
# 设置输出目录
output_folder = "I:\\TIF_res\\res_tif"
print("开始执行...")
arcpy.AddMessage("开始执行...")
# 获取所有Shapefile的名称（不包含扩展名）
shape_names = [os.path.splitext(os.path.basename(shp))[0] for shp in arcpy.ListFeatureClasses()]
# 遍历.tif文件
for raster_file in raster_files:
    town_name = os.path.splitext(raster_file)[0]  # 获取.tif文件的乡镇名称
    # 确保当前.tif对应的.shp文件存在
    if town_name in shape_names:
        # 构建.tif文件的完整路径
        raster_path = os.path.join(raster_folder, raster_file)
        # 构建输出栅格的完整路径
        out_raster_path = os.path.join(output_folder, f"Extract_{town_name}.tif")
        # 确保输出目录存在
        if not os.path.exists(output_folder):
            os.makedirs(output_folder)
        arcpy.AddMessage("正在执行处理" + town_name + "...")
        # 获取对应的.shp文件路径
        shp_file = os.path.join(arcpy.env.workspace, town_name + ".shp")
        # 在处理每个shapefile之前，设置环境的处理范围为该shapefile的范围
        arcpy.env.extent = arcpy.Describe(shp_file).extent
        arcpy.env.mask = shp_file
        try:
            if town_name != "":
                # 读取原始栅格的NoData值
                raster = arcpy.Raster(raster_path)
                original_nodata = raster.noDataValue
                # 提取掩膜区域
                out_extract_by_mask = arcpy.sa.ExtractByMask(raster, shp_file)
                # 保存提取后的栅格
                out_extract_by_mask.save(out_raster_path)
                # 将掩膜区域外的背景设为NoData
                if original_nodata is not None:
                    arcpy.management.SetRasterProperties(out_raster_path, nodata=f"1 {original_nodata}")
                arcpy.AddMessage(f"{town_name}按掩膜提取完成。")
        except Exception as e:
            arcpy.AddMessage(f"处理{town_name}时出错：{e}")
    else:
        arcpy.AddMessage(f"{town_name}对应的.shp文件未找到，跳过处理。")
print("所有匹配乡镇的按掩膜提取完成！")
arcpy.AddMessage("所有匹配乡镇的按掩膜提取完成！")

注意：

print()实在pyCharm控制台，arcpy.AddMessages是脚本信息输出

ArcGIS pro输出路径是地理数据库gdb的不需要文件后缀，输出路径是文件夹的，需要文件后缀