Python地理数据处理 十五：基于arcpy的批量操作

今天把最近遇到的一些问题和解决方法和大家进行分享，方便自己后面的回顾，也希望帮助大家解决自己的问题。后期如果大家需要的话，我会出一个关于Anusplin插值软件的一个详细教程，免费分享给大家。
arcpy可以在Arcgis自带的python2.7版本中使用，不需要再下载python。

1. 批量格式转换（grd to tiff）

# -*- coding: cp936 -*-
import arcpy
arcpy.env.workspace="E:\\ANUSSPLIN\\90DEM2\\"  #存放数据的文件夹
a=arcpy.ListRasters("*","grd")  #得到文件夹下所有的grd名称
for i in a:
    arcpy.RasterToOtherFormat_conversion(i,"E:\\ANUSSPLIN\\90DEM2\\tif\\","TIFF")

2. 批量定义投影(Albers)

import arcpy
arcpy.env.workspace = "E:\\ANUSSPLIN\\1000mTEM"
rasters = arcpy.ListRasters("*","tif")
Coordinate_System = "E:\\ANUSSPLIN\\1000mPRE\\coordinateSystem.tif"
for raster in rasters:
    print str(raster)
    arcpy.DefineProjection_management(raster,Coordinate_System)
print("OK！！！")

3. 批量投影转换(WGS 1984 to WGS 1984 Albers)

非常感谢大家的指正，我把这里再进行完善一下。

# -*- coding: cp936 -*-
import arcpy
arcpy.env.workspace = "E:/Work_GIS/practice" # 输入路径
rasters = arcpy.ListRasters("*","tif")
for raster in rasters:
    print str(raster)
    out = "E:/Work_GIS/practice/output/" + raster  # 输出路径
    arcpy.ProjectRaster_management(raster,out,"PROJCS['Asia_North_Albers_WGS84_LCR',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',105.0],PARAMETER['Standard_Parallel_1',25.0],PARAMETER['Standard_Parallel_2',47.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"
    )
print("ok!!!")

ProjectRaster_management中的投影信息在哪里可以找到：

结果显示：

3.1 转换前的栅格图像属性：

这里采用的重采样方法是最邻近法

3.2 转换后的栅格图像属性：

此处的采样方法默认为最邻近法

所以对栅格图像进行批量裁剪的时候，只需要改变投影坐标信息就行，重采样方法默认为最邻近法

4. ArcPy批量裁剪

# -*- coding: cp936 -*-
import arcpy
arcpy.CheckOutExtension("spatial")
arcpy.env.workspace = "E:\\ANUSSPLIN\\1000mPRE"#栅格文件路径
rasters = arcpy.ListRasters("*","tif")#也可以是其他的文件格式，如grd
mask = "E:\\China_map\\长株潭垃圾站点信息\\长沙.shp"
for raster in rasters:
    print str(raster)
    out = "E:\\ANUSSPLIN\\PRE_output\\" + raster
    arcpy.gp.ExtractByMask_sa(raster, mask, out)
    print("ma_" + raster + "has done!")
print("ok!!!")

4.1 常见错误

Traceback (most recent call last):
  File "D:\Datasets\arcgisDatasetFrost\clip\clip", line 20, in <module>
    print str(raster)
UnicodeEncodeError: 'ascii' codec can't encode characters in position 0-4: ordinal not in range(128)

   这个错误通常发生在尝试将非ASCII字符编码为ASCII字符集时，因为ASCII字符集仅支持128个字符，无法表示非ASCII字符。
  解决这个问题的一种方法是在打印或输出非ASCII字符之前，将其转换为正确的字符集，如UTF-8或Unicode。在第一行加入：

# -*- coding: utf-8 -*-

   但是有时候尝试了上面的方法之后还是不行，这里建议你将你要裁剪的栅格影像的文件名改为英文或数字（不要包含汉字），就可以解决99%的问题了。

5. ArcPy批量镶嵌

# -*- coding: cp936 -*-
import arcpy
import os
arcpy.CheckOutExtension("spatial")
base = r"E:\Download\2020年中国行政区划边界SHP数据-省、市\安徽省.shp"
#获取坐标系
out_coordSystem = arcpy.Describe(base).spatialReference
dataType = arcpy.Describe(base).DataType
piexlType = arcpy.Describe(base).pixelType
cellWidth = arcpy.Describe(base).meanCellWidth
bandCount = arcpy.Describe(base).bandCount
arcpy.env.workspace = "E:\\Download\\2020年中国行政区划边界SHP数据-省、市\\005-2020年中国行政区划边界-省、市-Shp\\2020年中国行政区划边界-省、市-Shp\\省边界（包含子区域）\\"
rasters = arcpy.ListRasters("*","shp")
outFolder = r"E:\Download\2020年中国行政区划边界SHP数据-省、市\output"
arcpy.MosaicToNewRaster_management(rasters, outFolder,"result.shp",
                                   out_coordSystem,"16_BIT_SIGNED",
                                   cellWidth,bandCount,"LAST","FIRST")