python之arcgis Query空间查询（点线面取交集）

皓空揽月

已于 2023-11-18 14:46:40 修改

阅读量176

点赞数

分类专栏： Python GIS 文章标签： arcgis python

于 2023-05-06 16:15:48 首次发布

本文链接：https://blog.csdn.net/weixin_38977885/article/details/130528126

版权

GIS 同时被 2 个专栏收录

2 篇文章 0 订阅

订阅专栏

Python

1 篇文章 0 订阅

订阅专栏

需求：通过arcgisserver根据绘制的几何图形区域(点、线、面)进行空间查询，返回的结果需要与绘制的几何图形区域进行裁剪取交集。

文章目录

1.python代码

#!/usr/bin/env python 
# -*- coding:utf-8 -*-
# 作者：zy
import json
import requests
import operator
from osgeo import ogr
from osgeo import osr


def arcgisQuery(request, folder, serviceName, layerId):
    # 判断是否是GET请求
    print(request.path)
    global spatialRel, where
    if request.method == 'GET':
        datas = request.GET
        where = datas.get('where')  # 属性查询条件
        geometry = datas.get('geometry')  # 空间查询集合
        geometryType = datas.get('geometryType')  # 空间类型
        spatialRel = datas.get('spatialRel')  # 空间关系
        outFields = datas.get('outFields')  # 输出字段
        returnGeometry = datas.get('returnGeometry')  # 返回几何图形
        returnCountOnly = datas.get('returnCountOnly')  # 返回总数
        pageNum = datas.get('pageNum')  # 页面返回大小
        pageSize = datas.get('pageSize')  # 页面返回大小
        orderByFields = datas.get('orderByFields')  # 排序字段


    # 初始化参数
    if pageNum is None:
        pageNum = 1
    if pageSize is None:
        pageSize = 10
    # mapserver空间关系转esri空间关系
    if spatialRel == 'Intersects':
        spatialRel = 'esriSpatialRelIntersects'
    elif spatialRel == 'Overlaps':
        spatialRel = 'esriSpatialRelOverlaps'
    elif spatialRel == 'Within':
        spatialRel = 'esriSpatialRelWithin'
    else:
        spatialRel = ''
	
	# geometry 入参这里是按照mapserver的传参格式，所以下面需要转换成arcgis支持的格式
    # mapserver几何图形转esri几何图形
    if geometryType == 'Point' and geometry is not None:
        geometry = geometry.replace(' ', ',')
        geometryType = 'esriGeometryPoint'
    elif geometryType == 'Polyline' and geometry is not None:
        split = geometry.split(',')
        outterlist = []
        innerlist = []
        for x in split:
            co = x.split()
            itlist = []
            itlist.append(float(co[0]))
            itlist.append(float(co[1]))
            innerlist.append(itlist)
        outterlist.append(innerlist)
        geometry = {
            "paths": outterlist
        }
        geometryType = 'esriGeometryPolyline'
    elif geometryType == 'Polygon' and geometry is not None:
        split = geometry.split(',')
        outterlist = []
        innerlist = []
        for x in split:
            co = x.split()
            itlist = []
            itlist.append(float(co[0]))
            itlist.append(float(co[1]))
            innerlist.append(itlist)
        outterlist.append(innerlist)
        geometry = {
            "rings": outterlist
        }
        geometryType = 'esriGeometryPolygon'
    else:
        pass

    # 设置url
    url = 'http://192.168.1.9:6080/arcgis/rest/services/' + folder + '/' + serviceName + '/MapServer/' + layerId + '/query'
    print(json.dumps(geometry))

    # 设置查询参数
    params = {
        "f": "json",
        "where": "1=1" if where is None else where,
        "geometryType": geometryType,
        "spatialRel": spatialRel,
        "geometry": json.dumps(geometry),
        "outFields": "*" if outFields is None else outFields,
        "returnGeometry": "true" if returnGeometry == "false" else returnGeometry
    }

    # 将GET请求发送到ArcGIS REST服务并获取响应
    response = requests.get(url, params=params)

    # 处理响应
    if response.status_code == 200 and 'error' not in response.json():
        count = 0
        total = 0
        newfeatures = []
        # 将响应解析为JSON
        data = response.json()
        # 从响应中获取要素
        features = data["features"]
        # print(data)
        mapLayer = {}
        # 创建列表以存储交叉点几何图形
        intersection_geometries = []
        # esri的几何图形转成geojson的几何图形
        geom = params['geometry']
        # 缓冲查询
        # if geometryType == 'esriGeometryPoint' or geometryType == 'esriGeometryPolyline':
        #     geom = geom.Buffer(0.00001)
        # 202305051559修改
        geom_of_geojson = esrijson_2_geojson(json.loads(geom), geometryType)
        wkid = data["spatialReference"]["wkid"]

        # 迭代要素并找到它们与矩形几何图形的交点
        for feature in features:
            fields = []
            intersectionArea = 0
            # 从要素中提取几何图形
            # geometry_json = json.dumps(feature["geometry"])
            geometry = feature["geometry"]
            # print('--', json.dumps(feature), '--')
            # 几何图形的esrijson转geojson
            geojson_geometry = esrijson_2_geojson(geometry, "esriGeometryPolygon")

            # 将 GeoJSON 字符串转换成 OGR Geometry 对象（几何对象），可以用于在 GIS 软件中进行空间分析、地理处理等操作
            geometry = ogr.CreateGeometryFromJson(str(geojson_geometry))

            # 定义feature
            newfeature = {"attributes": feature['attributes']}
            # 裁剪
            if geometry is not None:
                # 创建用于相交的矩形几何图形
                # rect_geometry_json = '{"type": "Polygon", "coordinates": [[[-123, 47], [-122, 47], [-122, 46], [-123, 46], [-123, 47]]]}'

                rect_geometry = ogr.CreateGeometryFromJson(json.dumps(geom_of_geojson))
                # 执行要素几何图形和矩形几何图形的相交
                # print(geometry)
                intersection_geometry = rect_geometry.Intersection(geometry)

                # 如果找到交叉点，请将其几何图形添加到交叉点几何图形阵列中
                if intersection_geometry is not None:
                    intersectionArea = get_intersection_area(intersection_geometry, wkid)
                    intersection_geometries.append(intersection_geometry.ExportToJson())
                    newfeature["geometry"] = intersection_geometry.ExportToJson()
                else:
                    intersectionArea = 0
                newfeature["intersectionArea"] = intersectionArea
                # 加入到features 中
                newfeatures.append(newfeature)
            else:
                print('rect_geometry is None')
        # 分页
        total = len(newfeatures)
        if pageNum is not None and pageSize is not None:
            newfeatures = pagination(int(pageNum), int(pageSize), newfeatures)

        # 封装图层对象
        mapLayer = {"fields": data['fields'], "features": newfeatures,
                    "geometryType": str.replace(geometryType, 'esriGeometry', ''), "total": total,
                    "pageNum": pageNum, "pageSize": pageSize}
        return mapLayer
    else:
        print("Error: request failed.")
        return response.json()


# esrijson转geojson的geometry
def esrijson_2_geojson(geometry, geometryType):
    if "esriGeometryPoint" == geometryType:
        coords = _get_coordinates1(geometry)
        geojson_geometry = {
            'type': 'Point',
            'coordinates': coords
        }
        return geojson_geometry
    elif "esriGeometryPolyline" == geometryType:
        coordinates = geometry['paths'][0]
        geojson_geometry = {
            'type': 'LineString',
            'coordinates': coordinates
        }
        return geojson_geometry
    elif "esriGeometryPolygon" == geometryType:
        coordinates = []
        for ring in geometry['rings']:
            ring_coords = _get_coordinates([ring])[0]
            # 将 "rings" 中的坐标顺序从逆时针方向转换为顺时针方向
            if _is_ccw(ring_coords):
                ring_coords.reverse()
            coordinates.append(ring_coords)
        geojson_geometry = {
            'type': 'Polygon',
            'coordinates': coordinates
        }
        return geojson_geometry


def _get_coordinates(paths):
    """将路径列表转换为坐标列表"""
    coordinates = []
    for path in paths:
        coord_list = [[x, y] for x, y in path]
        coordinates.append(coord_list)
    return coordinates

def _get_coordinates1(geometry):
    coordinates = []
    for s in geometry.split(','):
        coordinates.append(float(s))
    return coordinates

def _is_ccw(coords):
    """判断点集坐标顺序是否为逆时针方向"""
    area = sum(coords[i][0] * coords[i + 1][1] - coords[i + 1][0] * coords[i][1] for i in range(-1, len(coords) - 1))
    return area < 0


def pagination(pageNum, page_size, data_list):
    """
    :param pageNum: 当前页码
    :param page_size: 每页数据量
    :param data_list: 总数据列表
    :return: 当前页码对应的数据列表
    """
    start = (pageNum - 1) * page_size
    end = start + page_size
    return data_list[start:end]


def get_geometry_type(coordinates):
    """
    根据 coordinates 来判断几何对象的类型
    """
    if isinstance(coordinates, list) and len(coordinates) > 0:
        if all(isinstance(c, (float, int)) for c in coordinates):
            # 点类型
            return "Point"
        elif all(isinstance(p, list) and len(p) >= 2 and all(isinstance(c, (float, int)) for c in p) for p in
                 coordinates):
            # 线类型
            return "LineString"
        elif all(isinstance(ring, list) and len(ring) >= 4 and all(
                isinstance(p, list) and len(p) >= 2 and all(isinstance(c, (float, int)) for c in p) for p in ring) for
                 ring in coordinates):
            # 面类型
            return "Polygon"
    return None


# 交集面积
def get_intersection_area(intersection, eosg):
    # 使用 WGS 1984 空间参考计算面积
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromEPSG(eosg)
    intersection.AssignSpatialReference(spatial_ref)
    area = intersection.GetArea()
    # print('intersection--', float(area))
    return area


# 获取arcgis图例
def get_legend(request, floder, serviceName):
    # 获取参数
    if request.method == 'GET':
        # 设置请求 URL 和参数
        url = "http://192.168.1.9:6080/arcgis/rest/services/" + floder + "/" + serviceName + "/MapServer/legend"
        params = {
            "f": "json",  # 返回 JSON 格式
            "dpi": 96,  # 指定 DPI 像素密度
            "dynamicLayers": "",  # 默认情况下，将使用所有动态图层
        }
        # print(floder, serviceName)
        # 发送 GET 请求并处理响应
        response = requests.get(url, params=params)
        if response.status_code == 200:
            legend_json = response.json()  # 获取 JSON 数据
            # 处理JSON数据，例如获取图例符号、标签等信息
            # print(legend_json)
            return legend_json
        else:
            print("Failed to retrieve legend.")
            return None