之前,整理了 ArcGIS自定义坐标系统(以Albers为例),但这个只是改变显示的投影,而不是数据的投影。有时候需要转变数据的投影。
在ArcGIS中,可以使用投影栅格工具修改影像的投影。
如果需要批量投影的话,可以使用arcpy调用 ProjectRaster_management 函数。
首先构建投影文件,这里构建了-90~90°的Albers投影文件,具体参数是:
# False_Easting = 0.0
# False_Northing = 0.0
# Central_Meridian = 0.0
# Standard_Parallel_1 = i - 10.0
# Standard_Parallel_2 = i + 10.0
# Latitude_Of_Origin = 0.0
批量生成Albers投影文件:
for i in range(-90, 91, 1):
print i
# prjname = "Albers_0_50S"
prjtype = "Albers"
isNS = ""
if i > 0:
isNS = "N"
elif i < 0:
isNS = "S"
False_Easting = 0.0
False_Northing = 0.0
Central_Meridian = 0.0
Standard_Parallel_1 = i - 10.0
Standard_Parallel_2 = i + 10.0
Latitude_Of_Origin = 0.0
prjname = prjtype + "_" + str(Central_Meridian) + "_" + str(abs(i)) + isNS
prjstr = "PROJCS["
prjstr = prjstr + "\"" + prjname + "\","
prjstr = prjstr + "GEOGCS[\"GCS_WGS_1984\","
prjstr = prjstr + "DATUM[\"D_WGS_1984\","
prjstr = prjstr + "SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],"
prjstr = prjstr + "PRIMEM[\"Greenwich\",0.0],"
prjstr = prjstr + "UNIT[\"Degree\",0.0174532925199433]],"
prjstr = prjstr + "PROJECTION[\"" + prjtype + "\"],"
prjstr = prjstr + "PARAMETER[\"False_Easting\"," + str(False_Easting) + "],"
prjstr = prjstr + "PARAMETER[\"False_Northing\"," + str(False_Northing) + "],"
prjstr = prjstr + "PARAMETER[\"Central_Meridian\"," + str(Central_Meridian) + "],"
prjstr = prjstr + "PARAMETER[\"Standard_Parallel_1\"," + str(Standard_Parallel_1) + "],"
prjstr = prjstr + "PARAMETER[\"Standard_Parallel_2\"," + str(Standard_Parallel_2) + "],"
prjstr = prjstr + "PARAMETER[\"Latitude_Of_Origin\"," + str(Latitude_Of_Origin) + "],"
prjstr = prjstr + "UNIT[\"Meter\",1.0]]"
#print prjstr
prjfile = prjpath + prjname + ".prj"
fh = open(prjfile, 'w')
fh.write(prjstr)
fh.close()
如果想生成其他类型的投影文件,直接用文本打开.prj文件即可,修改部分参数就可以。
使用 ProjectRaster_management 函数修改投影,参数分别为:输入文件,输出文件,投影文件,重采样方式,和分辨率(米)
for i in range(-90, 91, 1):
inlltif = lltifpath + "lat_" + str(i) + ".tif"
if os.path.exists(inlltif) == False:
continue
prjname = prjtype + "_" + str(Central_Meridian) + "_" + str(abs(i)) + isNS
prjfile = prjpath + prjname + ".prj"
outaberls = tifAlberspath + "lat_" + str(i) + "_Aberls.tif"
arcpy.ProjectRaster_management(inlltif, outaberls, prjfile, "BILINEAR", "25000", "#", "#", "#")
投影前的影像(0.25°):
投影前的影像投影参数:
投影后的影像(25km):
投影后的影像投影参数:
欢迎大家批评指正!