以下代码是依据ChatGPT的答案，然后不断调整修改，并成功的代码案例。
-
需求：
这里有两个文件夹，分别存放10个栅格文件和3个矢量文件。 矢量文件与栅格文件都是地理坐标WGS84。首先根据矢量文件的空间位置确定对应的栅格文件。 其次，使用此矢量文件裁剪栅格文件。 第三，将裁剪后的栅格文件投影到 WGS84 UTM 中。 这里需要自动判断UTM条带号，是在南半球还是北半球。并且重采样到500m分辨率。第四，进行均值滤波。 以上请使用Python2.7中的Arcpy来完成。 -
代码：Vscode
import arcpy
import os
from arcpy.sa import *
import glob
# Set the workspace and get the list of raster and vector files
RasterFolder = "path/to/raster/folder"
raster_list = glob.glob(os.path.join(RasterFolder,"*.tif"))
arcpy.env.workspace = "path/to/vector/folder"
vector_list = arcpy.ListFeatureClasses()
# Set the output path for raster saving
Output_Folder = "path/to/raster/folder"
# Loop through the vector files
for vector_file in vector_list:
# Get the extent of the vector file
extent = arcpy.Describe(vector_file).extent
# Calculate the UTM zone
utm_zone = int((extent.XMin + 180) / 6) + 1
# Determine if the vector file is in the northern or southern hemisphere
if extent.YMin >= 0:
WKID = 32600 + utm_zone
utm_projection = arcpy.SpatialReference(WKID)
else:
WKID = 32700 + utm_zone
utm_projection = arcpy.SpatialReference(WKID)
# Loop through the raster files
for raster_file in raster_list:
# Get the extent of the raster file
raster_extent = arcpy.Describe(raster_file).extent
# Check if the extent of the vector file intersects with the raster file
if extent.disjoint(raster_extent):
continue
# clip: change into rectangle and clip img by rectangle
outFeature = "Rectangle_" + vector_file
arcpy.MinimumBoundingGeometry_management(vector_file,outFeature,"ENVELOPE","NONE")
out_extract = arcpy.sa.ExtractByMask(raster_file,outFeature)
out_cal = out_extract*(0.1)
# Project the clipped raster file into UTM projection and resample 500 by BILINEAR
output_projected_raster = vector_file[:-4] + "_projected.tif"
arcpy.ProjectRaster_management(out_cal, output_projected_raster, utm_projection,"BILINEAR", "500 500", "", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
#Spatial filter: Mean Focal Statistics
inRaster = arcpy.sa.FocalStatistics(output_projected_raster,"Rectangle 3 3 CELL","MEAN","DATA")
# Save result as Geotiff
outnamefocal = vector_file[:-4] + "_focal.tif"
outpathfocal = os.path.join(Output_Folder,outnamefocal)
inRaster.save(outpathfocal)
print(outpathfocal+"----------------Finish!")
# Clean up the temporary raster files
arcpy.Delete_management(outFeature)
arcpy.Delete_management(output_projected_raster)
- 体会：
-
整个代码过程中，耗费时间最多的是投影的步骤。这一步需要设置投影系统out_coor_system。可以从官网中知道，out_coor_system可以是这三类。
image.png
（1）坐标系的字符串表达形式（不推荐）
使用Arcmap中的ModelBuilder（拖拽工具，输入你的坐标，Export to Python Script）可以导出一个坐标系的字符串。北半球UTM的坐标字符串是如下表示：
"PROJCS['WGS_1984_UTM_Zone_" + str(utm_zone) + "N" + "',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',(6 * int(utm_zone[1:]) - 183)],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"
这里面需要改动的就是条带号、南北半球、'False_Northing'、中央经度。但是这个方法运行，一直在报错。
（2）prj文件（推荐）
首先，必须要有这么一个prj文件，才可以使用。如果只需要投影一个地方，比如200张都是投影UTM 50N， 可以找一个UTM 50N或者自己创建一个，再进行投影。但是如果要投影多个地方，那么这个方法不推荐。因为要一个一个建立prj文件，也是挺累。
（3）SpatialReference对象（推荐）
关于这个可以看这篇微信文章：
https://mp.weixin.qq.com/s/nURoQhZ808YkhSc0Ep68bg
Arcmap所有坐标的WKID与名称：
https://desktop.arcgis.com/zh-cn/arcmap/latest/analyze/arcpy-classes/pdf/projected_coordinate_systems.pdf
至于UTM的WKID，我是根据Arcmap的文件，推测出以下这样的规律：
# Determine if the vector file is in the northern or southern hemisphere
if extent.YMin >= 0:
WKID = 32600 + utm_zone # North
utm_projection = arcpy.SpatialReference(WKID)
else:
WKID = 32700 + utm_zone # South
utm_projection = arcpy.SpatialReference(WKID)
结果跑出来是对应上的。
网友评论