arcpy获取国土用地坐标范围的几种尝试

最新推荐文章于 2024-05-31 09:00:00 发布

规划酱

最新推荐文章于 2024-05-31 09:00:00 发布

阅读量2k

点赞数 1

分类专栏：实用技巧文章标签： python arcgis

本文链接：https://blog.csdn.net/fenfee/article/details/108784270

版权

实用技巧专栏收录该内容

25 篇文章 10 订阅

订阅专栏

在做国土空间规划的时候，准备在GIS中爬取对应范围的poi点，首先得获取整个用地文件的坐标范围。
最简单的办法是看图层属性，图层属性>源>范围，在GCGCS2000地里坐标系里，范围就是经纬度坐标。但要自动获取问题就来了！
在这里插入图片描述

1、获取图层框的范围

由于在用地图中，采用的是投影坐标，整个图层框默认就成了投影坐标，
图层框的范围获取用extent属性，出来的结果是大数的投影坐标，单位是m

mxd = arcpy.mapping.MapDocument("CURRENT") #获取当前map
df  = arcpy.mapping.ListDataFrames(mxd)[0] #获取第一个DataFrame
print(df.extent)
>>>40432311.6013862 3537106.75526 40440489.6550138 3543035.84414 NaN NaN NaN NaN

2、获取图层的范围

接着尝试用获取图层的范围，用layer.getextent()属性，特地用的地里坐标

ly_name = '三调' #用地里坐标系的图层
mxd = arcpy.mapping.MapDocument("CURRENT")
layer =  arcpy.mapping.ListLayers(mxd, ly_name)[0]
print(layer.getExtent(False)) #用False才是获取范围
>>>40434191.8945 3537376.2593 40438609.3619 3542766.3401 NaN NaN NaN NaN

这里出来的还是投影坐标，是地里坐标被转换成投影做的范围.

3、用游标获取质心坐标

GIS提供了获取几何要素属性的游标，帮助里是这么说的：

几何属性可通过在字段列表中指定令牌 SHAPE@ 进行访问。
使用 For 循环可迭代搜索游标。搜索游标也支持 With 语句；使用 With 语句可保证数据库锁的关闭和释放，并重置迭代。
利用属性条件和/或空间条件，可以限制 SearchCursor 返回的记录。
使用 SHAPE@ 访问整个几何是一种代价较高的操作。如果仅需要几何信息（如点的 x、y 坐标），请使用 SHAPE@XY、SHAPE@Z 和 SHAPE@M 之类的令牌进行更快速、高效的访问。
以令牌（如 OID@）取代字段名称可访问更多的信息：
SHAPE@XY —一组要素的质心 x,y 坐标。
SHAPE@TRUECENTROID —一组要素的真正质心 x，y 坐标。
SHAPE@X —要素的双精度 x 坐标。
SHAPE@Y —要素的双精度 y 坐标。
SHAPE@Z —要素的双精度 z 坐标。
SHAPE@M —要素的双精度 m 值。
SHAPE@JSON — 表示几何的 esri JSON 字符串。
SHAPE@WKB —OGC 几何的熟知二进制 (WKB) 制图表达。该存储类型将几何值表示为不间断的字节流形式。
SHAPE@WKT —OGC 几何的熟知文本 (WKT) 制图表达。其将几何值表示为文本字符串。
SHAPE@ —要素的几何对象。
SHAPE@AREA —要素的双精度面积。
SHAPE@LENGTH —要素的双精度长度。
OID@ —ObjectID 字段的值。

#质心
import time
s_time=time.time()
ly_name = '三调' #用地里坐标系的图层
list_x = []
list_y = []
with arcpy.da.SearchCursor(ly_name,["SHAPE@XY"]) as cursor:
    for row in cursor:
        list_x.append(row[0][0])
        list_y.append(row[0][1])
list_x.sort() #把所有质心X坐标排个序，最大最小就出来了
list_y.sort()
xmin = list_x[0]
xmax = list_x[-1]
ymin = list_y[0]
ymax = list_y[-1]
e_time=time.time()
print("程序共耗时： %.2f s!" % (e_time - s_time))
print(ymin,xmin,ymax,xmax)
>>>程序共耗时： 0.08 s!
>>>(31.95860735036497, 119.30386817778967, 32.006091115992284, 119.35003990813219)

因为基于三调的用地是几千个图斑组成的，用质心的XY来近似的获取了范围。
再写个代码看看大小对不对

#画个范围框
def kuang(ymin,xmin,ymax,xmax,out_path):
    kuang = [[xmin,ymin],[xmin,ymax],[xmax,ymax],[xmax,ymin]] #成为几个点
    coordinate = arcpy.SpatialReference(4490) #GCGCS2000坐标
    arcpy.CreateFeatureclass_management (out_path, 'kuang',"POINT", "", "","", coordinate) #新建个点图层
    rows=arcpy.InsertCursor('kuang')
    for i in range(len(kuang)):
        row = rows.newRow() #新建行
        vertex = arcpy.CreateObject("Point") #新建点
        vertex.X = kuang[i][0]
        vertex.Y = kuang[i][1]
        row.shape = vertex
        rows.insertRow(row)    #插入行
    arcpy.PointsToLine_management("kuang", out_path+'\\lines',"","" ,"CLOSE" ) #点要素连线
out_path = layer.workspacePath
kuang(31.95860735036497, 119.30386817778967, 32.006091115992284, 119.35003990813219,out_path)

看看框
在这里插入图片描述
范围是不全的，切了小部分地块图斑。

4、用游标获取几何边界点坐标

还是那个游标，获取所有几何元素，遍历元素的部分，然后遍历每个点，如果存在，获取点坐标，获取第一个点之后开始比大小，最终获得最大和最小的，笨办法，理论上肯定能出结果

#所有几何要素边界点比大小
s_time=time.time()
n = 1
with arcpy.da.SearchCursor(output, ["OID@", "SHAPE@"]) as cursor:
    for row in cursor:
    # 分解每个要素的部分（比如环，有内部和外部之分）
        for part in row[1]:
        #每个点
            for pnt in part:
                if pnt:
                    if n == 1:
                        xmin = pnt.X
                        xmax = pnt.X
                        ymin = pnt.Y
                        ymax = pnt.Y
                        n = 0
                    if pnt.X < xmin : xmin = pnt.X
                    if pnt.X > xmax : xmax = pnt.X
                    if pnt.Y < ymin : ymin = pnt.Y
                    if pnt.Y > ymax : ymax = pnt.Y
e_time=time.time()
print("程序共耗时： %.2f s!" % (e_time - s_time))
print(ymin,xmin,ymax,xmax)
>>>程序共耗时： 5.73 s!
>>>(31.957932880000044, 119.30361648600001, 32.006381426000075, 119.35055022000006)

画个框看下，这下全包进去了，就是速度比较慢，改良一下
在这里插入图片描述

#所有点排序
s_time=time.time()
list_x = []
list_y = []
with arcpy.da.SearchCursor(output, ["OID@", "SHAPE@"]) as cursor:
    for row in cursor:
    # 分解每个要素的部分（比如环，有内部和外部之分）
        for part in row[1]:
        #每个点
            for pnt in part:
                if pnt:
                    list_x.append(pnt.X)
                    list_y.append(pnt.Y)
list_x.sort()
list_y.sort()
xmin = list_x[0]
xmax = list_x[-1]
ymin = list_y[0]
ymax = list_y[-1]
e_time=time.time()
print("程序共耗时： %.2f s!" % (e_time - s_time)) 
print(ymin,xmin,ymax,xmax)
>>>程序共耗时： 4.31 s!
>>>(31.957932880000044, 119.30361648600001, 32.006381426000075, 119.35055022000006)

稍微快了一点，还是不理想。

5、用游标获取几何要素范围

# 几何要素范围
s_time=time.time()
list_x = []
list_y = []
with arcpy.da.SearchCursor(output, 'SHAPE@') as cursor:
    for row in cursor:
        extent = row[0].extent
        extent = str(extent).split()
        list_x.append(extent[0])  
        list_x.append(extent[2]) 
        list_y.append(extent[1]) 
        list_y.append(extent[3])
list_x.sort()
list_y.sort()
xmin = list_x[0]
xmax = list_x[-1]
ymin = list_y[0]
ymax = list_y[-1]
e_time=time.time()
print("程序共耗时： %.2f s!" % (e_time - s_time)) 
print(ymin,xmin,ymax,xmax)
>>>程序共耗时： 0.29 s!
>>>('31.95793288', '119.303616486', '32.006381426', '119.35055022')

这个速度还是比较理想的

6、方法5比较理想，完整测试一下

最后包装一下，两个思路获取shp坐标范围，
一是把用地的投影坐标通过投影工具转地里坐标，然后获取多几何元素的shp范围坐标
二是把用地用融合工具融合成一个几何元素，再把投影坐标通过投影工具转地里坐标，然后单几个元素获取shp范围

# 单几何元素的shp范围
def get_extent(layer):
    '''
    layer 图层名或者是路径
    '''
    with arcpy.da.SearchCursor(layer, 'SHAPE@') as cursor:
        for row in cursor:
            extent = row[0].extent
        ex = str(extent).split()         
    return ex[1]+','+ex[0]+','+ex[3]+','+ex[2]

# 多几何元素的shp范围
def get_extents(layer):
    '''
    layer 图层名或者是路径
    '''
    list_x = []
    list_y = []
    with arcpy.da.SearchCursor(layer, 'SHAPE@') as cursor:
        for row in cursor:
            extent = row[0].extent
            extent = str(extent).split()
            list_x.append(extent[0])  
            list_x.append(extent[2]) 
            list_y.append(extent[1]) 
            list_y.append(extent[3])
    list_x.sort()
    list_y.sort()
    return list_y[0]+','+list_x[0]+','+list_y[-1]+','+list_x[-1]

def outpath(): #存储路径
    mxd = arcpy.mapping.MapDocument("CURRENT")
    sd_layer =  arcpy.mapping.ListLayers(mxd, ly_name)[0]
    path = sd_layer.workspacePath
    if path[-3:] == u'gdb'or path[-3:] == u'mdb':
        return path+'\\rong',path+'\\dili'        
    else:
        return path+'\\rong.shp',path+'\\dili.shp'

#融合并将投影坐标转GCG2000地理坐标
def ds_project(sc_code = 4490): 
    '''
    转换坐标系代码
    4326 GCS_WGS_1984 
    4490 GCS_China_Geodetic_Coordinate_System_2000
    4528 CGCS2000_3_Degree_GK_Zone_40
    '''
    coordinate  = arcpy.SpatialReference(sc_code) 
    arcpy.Dissolve_management(ly_name, outpath1) #融合
    arcpy.Project_management(outpath1, outpath2, coordinate) #投影

#将投影坐标转GCG2000地理坐标
def project(sc_code = 4490): 
    '''
    转换坐标系代码
    4326 GCS_WGS_1984 
    4490 GCS_China_Geodetic_Coordinate_System_2000
    4528 CGCS2000_3_Degree_GK_Zone_40
    '''
    coordinate  = arcpy.SpatialReference(sc_code) 
    arcpy.Project_management(ly_name, outpath2, coordinate) #投影

测试一下：

#思路一：投影转坐标，获取范围
if __name__ == '__main__':
	s_time=time.time()
	ly_name = '三调'
	#路径
	outpath1,outpath2 = outpath()
	#融合并转 CS_China_Geodetic_Coordinate_System_2000
	project()
	#多几何元素范围
	extent = get_extents('dili')
	print(extent)
	e_time=time.time()
	print("耗时： %.2f s!" % (e_time - s_time))
>>>31.95793288,119.303616486,32.006381426,119.35055022
>>>耗时： 6.27 s!

#思路二：融合，投影转坐标，获取范围
if __name__ == '__main__':
	s_time=time.time()
	ly_name = '三调'
	#路径
	outpath1,outpath2 = outpath()
	#融合并转 CS_China_Geodetic_Coordinate_System_2000
	ds_project()
	#单几何元素范围
	extent = get_extent('dili')
	print(extent)
	e_time=time.time()
	print("耗时： %.2f s!" % (e_time - s_time))
>>>31.95793288,119.303616486,32.006381426,119.35055022
>>>耗时： 8.36 s!

看来都不是太理想啊，还是挺慢的！
有更好的办法请留言，我继续尝试别的方法。

规划酱

关注

1
点赞
踩
7

收藏

觉得还不错? 一键收藏
打赏
2
评论
arcpy获取国土用地坐标范围的几种尝试

在做国土空间规划的时候，准备在GIS中爬取对应范围的poi点，首先得获取整个用地文件的坐标范围。最简单的办法是看图层属性，图层属性>源>范围，在GCGCS2000地里坐标系里，范围就是经纬度坐标。但要自动获取问题就来了！1、获取图层框的范围由于在用地图中，采用的是投影坐标，整个图层框默认就成了投影坐标，图层框的范围获取用extent属性，出来的结果是大数的投影坐标，单位是mmxd = arcpy.mapping.MapDocument("CURRENT") #获取当前mapdf
复制链接

扫一扫