Arcpy常用代码
1.设置环境变量
# -*- coding: utf-8 -*-：
import arcpy
envfile = r"D:\Work\doing\CNN_RNN\data\beijing\inputdata\ALT.tif"
prjfile = r"D:\Work\doing\CNN_RNN\data\beijing\origindata\boundary.prj"
arcpy.env.overwriteOutput = 1 # 可覆盖
arcpy.env.workspace = r"D:\Work\doing\CNN_RNN\data\processdata" # 设置工作路径
arcpy.env.outputCoordinateSystem = prjfile
arcpy.env.extent = "MAXOF" # 设置范围
# arcpy.env.extent = arcpy.Extent(454871.698900, 4253992.158900, 624251.698900, 4437472.158900) # "XMin, YMin, XMax, YMax"
arcpy.env.extent = envfile
# arcpy.env.cellSize = 30 # 像元大小
arcpy.env.cellSize = envfile
arcpy.env.snapRaster = envfile
2. 扇形等分圆
# -*- coding: utf-8 -*-
import arcpy
import os
arcpy.env.overwriteOutput = 1 # 可覆盖
# --------------- 自定义参数 -------------------
# 圆的等分个数
n = 8
# 中心点坐标
x = 534874.483
y = 4309287.494
max_dis = 80
step = 5
# 圆形半径的距离，单位：km
distances = [j for j in range(step,max_dis,step)] # 135最大距离，5间隔，单位km
# distances = [80]
# 工作目录
env = arcpy.env.workspace = r"D:\Work\doing\00urbanAgglomerations\Result\data\beijing" + "/"
# tif遥感图路径
tif = r"D:\Work\doing\00urbanAgglomerations\figure\English\09Discuss\data\deeplearn\p1_plan1_beijing.tif"
output = "sector_{}_{}.shp".format(str(step),str(max_dis)) # 切割后最终结果
outFeatureClass = "circle_{}_{}.shp".format(str(step),str(max_dis)) # 切割前
# ---------------------------------------------
# 读取tif空间参考信息
tifsr = arcpy.Describe(tif).spatialReference
arcpy.env.outputCoordinateSystem = tifsr
fc = "point.shp"
if os.path.exists(env + fc):
arcpy.Delete_management(fc)
arcpy.CreateFeatureclass_management(env,fc,'POINT')
arcpy.AddField_management(fc, "NAME", "TEXT", "", "", "", "", "", "", "")
xy_values = [('point1', (x, y))]
cursor = arcpy.da.InsertCursor(fc, ("NAME", "SHAPE@XY"))
for row in xy_values:
cursor.insertRow(row)
del cursor
inFeatures = fc
bufferUnit = "Kilometers"
if os.path.exists(env + outFeatureClass):
arcpy.Delete_management(outFeatureClass)
arcpy.MultipleRingBuffer_analysis(inFeatures, outFeatureClass, distances, bufferUnit, "", "ALL")
degrees = [(j + 1) * 360 / n - (180/n) for j in range(n)]
if os.path.exists(env + "table.dbf"):
arcpy.Delete_management("table.dbf")
tablefc = arcpy.CreateTable_management(env, "table.dbf")
arcpy.AddField_management(tablefc, "distance", "DOUBLE", "", "", "", "", "", "", "")
arcpy.AddField_management(tablefc, "x", "DOUBLE", "", "", "", "", "", "", "")
arcpy.AddField_management(tablefc, "y", "DOUBLE", "", "", "", "", "", "", "")
arcpy.AddField_management(tablefc, "degrees", "FLOAT", "", "", "", "", "", "", "")
cursor = arcpy.InsertCursor(tablefc)
i = 0
while i < len(degrees):
row = cursor.newRow()
row.distance = max_dis + 1
row.x = x
row.y = y
row.degrees = degrees[i]
cursor.insertRow(row)
i += 1
del cursor
arcpy.BearingDistanceToLine_management("table.dbf", "lines", "x", "y", "distance", bufferUnit, "degrees",
"DEGREES", "", "", tifsr)
if os.path.exists(env + output):
arcpy.Delete_management(output)
arcpy.FeatureToPolygon_management(["lines.shp", outFeatureClass], output)
# ----------------删除中间数据------------------
# if os.path.exists(env + outFeatureClass):
# arcpy.Delete_management(outFeatureClass)
if os.path.exists(env + "lines.shp"):
arcpy.Delete_management("lines.shp")
if os.path.exists(env + "info"):
arcpy.Delete_management("info")
# if os.path.exists(env + "point.shp"):
# arcpy.Delete_management("point.shp")
if os.path.exists(env + "table.dbf"):
arcpy.Delete_management("table.dbf")
# ---------------------------------------------
不用代码的方法：http://www.360doc.com/content/12/0408/11/6743016_201876186.shtml
3. 分区统计导出excel
import arcpy
from arcpy import env
from arcpy.sa import *
import os
arcpy.CheckOutExtension("spatial")
# Set environment settings
env.workspace = r"G:\01UrbanWaterSourceSecurity" # 路径
arcpy.env.overwriteOutput = 1
# Set local variables
inZoneData = "data/shp/basin.shp" # shp
zoneField = "ID" # 字段
names = ['Chen', 'Harmonized_DN_NTL_2020_simVIIRS', 'LUH_crops_m_fertl_2020'] # tif
if not os.path.exists("result/dbf/"):
os.makedirs("result/dbf/")
if not os.path.exists("result/excel/"):
os.makedirs("result/excel/")
for name in names:
inValueRaster = "data/tif/" + name + ".tif"
outTable = "result/dbf/" + name + ".dbf"
out_xls = "result/excel/" + name + ".xls"
# Execute ZonalStatisticsAsTable
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster,
outTable, "DATA", "SUM")
arcpy.TableToExcel_conversion(outTable, out_xls)
print(inValueRaster)
4. 裁剪(做到结果行列行一致)
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")
oundary = r"D:\Work\doing\00urbanAgglomerations\data\origindata\boundary.shp" # 范围shp
envfile = r"D:\Work\doing\00urbanAgglomerations\data\inputdata\ALT.tif" # 标准数据
prjfile = r"D:\Work\doing\00urbanAgglomerations\data\origindata\boundary.prj"# 投影
processpath = r"D:\Work\doing\00urbanAgglomerations\data\processdata"# 中间过程数据存放路径
arcpy.env.overwriteOutput = 1 # 可覆盖
arcpy.env.workspace = processpath # 设置工作路径
arcpy.env.outputCoordinateSystem = prjfile
arcpy.env.extent = "MAXOF" # 设置范围
# arcpy.env.extent = arcpy.Extent(454871.698900, 4253992.158900, 624251.698900, 4437472.158900) # "XMin, YMin, XMax, YMax"
arcpy.env.extent = envfile
arcpy.env.cellSize = 500 # 像元大小
arcpy.env.cellSize = envfile
arcpy.env.snapRaster = envfile
# Local variables:
Parallel_Processing_Factor = "0"
inMaskData = r"D:\Work\doing\00urbanAgglomerations\data\origindata\envelop.shp"
path1 = r'D:\Work\data\CLCD'
path2 = r"D:\Work\doing\00urbanAgglomerations\data\origindata\lulc_clcd"
names = [1992,1993,1995,1996,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2012,2015,2020]
for name in names:
inputdir = path1 + "/" + "CLCD_v01_" + str(name) + ".tif"
outputdir = path2 + "/" + "lulc" + str(name) + ".tif"
outExtractByMask = ExtractByMask(inputdir, inMaskData)
outExtractByMask.save(outputdir)
print(name)
print("OK")
4. 添加字段并按属性赋值
# -*- coding:utf-8 -*-
import arcpy
import sys,locale
reload(sys)
sys.setdefaultencoding("utf-8")
arcpy.env.overwriteOutput = 1
Parallel_Processing_Factor = "0"
path = r"D:\Work\doing\00urbanAgglomerations\Result\data"
arcpy.env.workspace = path
def getNewValue(value):
if value == 0:
i = 1
elif value > 0 and value <= 50:
i = 2
elif value > 50:
i = 3
return i
#设置目标数据
path = r"D:\Work\doing\00urbanAgglomerations\Result\data"
names_new = ['urban2020', 'ssp5', 'flus_beijing_2020', 'flus_shijiazhuang_2020',]
for i in range(0, len(names_new)):
data = path + "/LEI/new/" + names_new[i] + ".shp"
print(data)
#添加字段type
arcpy.AddField_management(data,'type','SHORT')
#根据LEi字段逐行给type赋值
rows = arcpy.UpdateCursor(data)
for row in rows:
row.type = getNewValue(row.LEI)
rows.updateRow(row)
5. 按时shp属性裁剪tif
# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
arcpy.env.overwriteOutput = 1 # 可覆盖
shp =r"D:\Work\doing\00urbanAgglomerations\Result\data\direction\boundary.shp"#要处理的shp文件路径
out_path=r"D:\Work\doing\00urbanAgglomerations\Result\data\direction\ssp5"#输出文件夹的路径
input_tif = r"D:\Work\doing\00urbanAgglomerations\Result\data\direction\sim\ssp5.tif"
FieldsName = "name" # 选择的裁剪字段
arcpy.env.workspace = out_path
with arcpy.da.SearchCursor(shp, ["SHAPE@",FieldsName]) as cursor:
for row in cursor:
shp_name=str(row[1])+'.shp'#输出文件名 确保均为字符串类型，注意命名是唯一的，这样就能将该shp数据中的所有要素都导出
arcpy.FeatureClassToFeatureClass_conversion (row[0],out_path,shp_name)
shp_path = out_path + "/" + shp_name
outputdir = out_path + "/" + str(row[1]) + ".tif"
outExtractByMask = ExtractByMask(input_tif, shp_path)
outExtractByMask.save(outputdir)
arcpy.Delete_management(shp_path)
print(outputdir)
浙公网安备 33010602011771号