这个可以用arcpy做,但太麻烦了,想了个简单的法子,直接用属性计算器计算了。
我用的是10.2版本的arcgis,而且是没汉化的英文版
主要分两部分,一部分是提取图斑包络四边形顶点,另一个是计算图幅号
转地理坐标系
如果发现计算的坐标不是经纬度,那得转地理坐标系。我转的是CGCS2000,大家按照数据本来的坐标系转对应的坐标系
不想转的话也可以,直接改图层frame属性
到时候计算的时候选下面的选项:
提取图斑包络四边形顶点
方法参考于:liuxupiaoshi的博客文章
第一步、新建一个属性ID 直接从零开始赋值(等于FID就可以了)
第二步、生成包络图形
第三步、提取特征点
第四步、新建X 、 Y 属性 ,计算点坐标
注意计算的时候坐标系别搞错了
第五步、把点特征叠加到包络矩形,并作汇总
这里用的是相交工具:
对结果用提前设置好的ID作汇总
输出的是一个dbf表格
第六步、把这个表格用接回去
在原始图斑的属性表内利用ID属性进行join操作
结果如下
计算图幅号
其中,图幅号计算方法思路参考了darlun的博客文章中的代码
第一步、新建多个属性字段:
除了那个计数的属性,其他都是TEXT属性的
第二步,计算四个角点的图幅号,并用逗号隔开输出到 “总输出“ 属性
打开Field Calculator
开始写脚本
LatCharArray = ['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V']#图幅比例尺 scaleArray = [1000000,500000,250000,100000,50000,25000,10000,5000]#图幅纬差,单位秒 latDArray = [14400,7200,3600,1200,600,300,150,75]#图幅经差,单位秒 lonDArray = [21600,10800,5400,1800,900,450,225,112.5]SCALE_1000000 = 0SCALE_500000 = 1SCALE_250000 = 2SCALE_100000 = 3SCALE_50000 = 4SCALE_25000 = 5SCALE_10000 = 6SCALE_5000 = 7#lat 为 纬度, lon为经度 1:5000的scaleID为7def getSheetNumber(lat,lon,scaleID):f = lat*3600;#以秒表示的纬度r = lon*3600;#以秒表示的经度numString = ''a = int(f/(4*3600)); #纬度序号b = int(r/(6*3600))+31; #经度序号 c = int(4*3600/latDArray[scaleID])-int((f%(4*3600))/latDArray[scaleID]); #在1:1000000万图幅中的序号 纬度 d = int((r%(6*3600))/lonDArray[scaleID])+1; #在1:1000000万图幅中的序号 经度numString += LatCharArray[a];numString += ""+ str(b);#1:1000000万if(scaleID==0):return numString; numString +=LatCharArray[scaleID];if c<10: numString += "00"numString += str(c)elif c<100: numString += "0"numString += str(c)elif c<1000:numString += str(c) else: return "" if d<10: numString += "00"numString += str(d) elif d<100:numString += "0" numString += str(d) elif d<1000: numString += str(d); else:return ""return numStringdef f(xmin,xmax,ymin,ymax):return getSheetNumber(ymax,xmax,7)+ ',' + getSheetNumber(ymax,xmin,7)+ ',' + getSheetNumber(ymin,xmax,7)+ ',' + getSheetNumber(ymin,xmin,7)# F1 F2 是我自己写的小代码,# 补零位def F1(a):if a<10:strout = '00' + str(int(a))elif a < 100:strout = '0' + str(int(a))else:strout = str(int(a))return stroutdef F2(x,y):#区域编号头部str1 = "I50H"#1:5000标准经差d1 = 1.0/60 + 52.5/3600#标准纬度差d2 = 1.0/60 + 15.0/3600c = 4/d2 - int(( y % 4) /d2)d = int((x % 6) / d1) + 1str2 = F1(c)str3 = F1(d)str4 = str1 + str2 + str3return str4
用field计算器开始运行代码,把代码粘到对应的框子里面
结果
第三步,计算重复与输出结果的中间步骤
中间输出格式为 [跨图幅数],[所有的图幅号]
代码是
#中间def f(str0):str0 = str0.split(',')a = str0[0]b = str0[1]c = str0[2]d = str0[3]count = 1str1 = aif b!=a:count +=1str1 =str1 + ',' + bif c !=a and c != b:count +=1str1 =str1 + ',' + cif d != a and d != b and d != c:count +=1str1 =str1 + ',' + dreturn str(count) + ',' + str1
第四步,输出结果
分别提取跨图幅数和图幅号
#countdef f (a):return int(a[0])
#图幅号def f (a):return a[2:]
结果如下
第五步,修饰
解除join
按照不同的图幅号显示颜色;欧克 完成了