python 地理分布图 气泡图 python地理空间分析

5.6.1点包容性公式

主要利用光影投射法执行检查操作。
该方法会从测试点创建一条直线并穿过多边形，之后会计算其和多边形每条边相交后产生的点的个数。

• 如果该数目是偶数，那么点在多边形外部；
• 如果该数目是奇数，那么点在多边形内部。

def point_in_poly(x, y, poly):
    # 检查是不是顶点
    if (x, y) in poly:
        print("是顶点")
        return True


    # 检查影像是否在边界
    length = len(poly)
    for i in range(length):
        p1 = None
        p2 = None
        if i == 0:
            p1 = poly[0]
            p2 = poly[1]
        else:
            p1 = poly[i-1]
            p2 = poly[i]
        # 判断点在poly范围内
        if p1[1] == p2[1] and p1[1] == y and x > min(p1[0],p2[0]) and x < max(p1[0], p2[0]):
            print("在边界上")
            return True

    n = len(poly)
    inside = False

    p1x, p1y = poly[0]
    for i in range(n+1):
        p2x, p2y = poly[i % n]
        #test = i % n
        #print("test:",test)
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x, p1y = p2x, p2y

    if inside:
        print("在内部")
        return True
    print("在外部")
    return False

if __name__ == '__main__':

    #  判断一个点是否包含在某区域
    myPolygon = [(-70.593016, -33.416032), (-70.589604, -33.415370),
                 (-70.589046, -33.417340), (-70.592351, -33.417949),
                 (-70.593016, -33.416032)]
    # 测试位置点1
    lon = -70.592000
    lat = -33.416000

    print(point_in_poly(lon, lat, myPolygon))


    # 测试位置点2
    lon = -70.593016
    lat = -33.416032
    print(point_in_poly(lon, lat, myPolygon))
    # 测试位置点3
    lon = -73.593016
    lat = -33.416032
    print(point_in_poly(lon, lat, myPolygon))

测试位置点1：
在内部
True

测试位置点2：
是顶点
True

测试位置点3：
在外部
False

5.6.2边框查询

使用个简单的边框将一个复杂的特征集子集化，之后将其保存为一个新的Shapefile文件。

"""Shapefile spatial query"""

# https://github.com/GeospatialPython/Learn/raw/master/roads.zip

import shapefile

# 读取数据
r = shapefile.Reader(r"roads\roadtrl020")
# 创建一个w文件，用于写入在选择区内的r数据
with shapefile.Writer(r"roads\Puerto_Rico_Roads", r.shapeType) as w:
    # 字段
    w.fields = list(r.fields)
    # 定义选择区域的范围
    xmin = -67.5
    xmax = -65.0
    ymin = 17.8
    ymax = 18.6
    # 循环读取r文件的记录
    for road in r.iterShapeRecords():
        # 定义shape
        geom = road.shape
        # D定义记录
        rec = road.record
        # 得到每条记录的外包矩形范围
        sxmin, symin, sxmax, symax = geom.bbox
        # 判断是否在所选的范围内
        if sxmin < xmin:
            continue
        elif sxmax > xmax:
            continue
        elif symin < ymin:
            continue
        elif symax > ymax:
            continue
        # 在范围内的线被添加到新的图层中
        w.line([geom.points])
        w.record(*list(rec))

roadtrl020数据

右下角是要查询的结果

5.6.3 属性查询一

使用属性表获取矢量数据子集的方法，通过记录表中的城市人口密度来进行筛选，选择城市人口密度小于5000的要素。

import shapefile
# 读取矢量文件
r = shapefile.Reader(r"MS_UrbanAnC10/MS_UrbanAnC10.shp")
with shapefile.Writer(r"MS_UrbanAnC10/MS_UrbanAnC10_subset.shp", r.shapeType) as w:
    w.fields = list(r.fields)

    selection = []
    # 打印第一条记录
    print(r.record())
    # 打印所有记录
    print(r.records())
    # 注意enumerate的用法，
    for rec in enumerate(r.records()):
        if rec[1][15]<5000:
            # print(rec[0])
            # print(rec[1])
            # print(rec[1][15])
            selection.append(rec)

    for rec in selection:
        # 添加形状
        w.poly([r.shape(rec[0]).points])
        # 添加记录
        w.record(*list(rec[1]))

结果：

5.6.3 属性查询二

使用fiona库在上一个示例中的表现，该库能够方便地调用OGR库。

import fiona

with fiona.open(r"MS_UrbanAnC10/MS_UrbanAnC10.shp") as sf:
    #将使用嵌套语句适当缩减打开和关闭文件代码的数量。
    filtered = filter(lambda f:f["properties"]["POP"]<5000,sf)
    #shapefile文件格式驱动
    drv = sf.driver
    # 参考坐标系
    crs = sf.crs
    #Dbf架构
    schm = sf.schema
    #定义子集文件名
    subset = r"MS_UrbanAnC10/MS_UrbanAnC10_FionaSubset.shp"
    with fiona.open(subset,"w",driver = drv,crs = crs,schema = schm) as w:
        for rec in filtered:
            w.write(rec)

---

总结

本节介绍了两种空间查询，两种属性查询的案例。