# Author：hhx
# Date：2022/8/12 10:14
# Description：在多个面属性的shp文件里面选点，判断是否在shp内


import shapefile
from shapely import geometry
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib
import numpy as np
import pylab as plt
from osgeo import osr
from tqdm import tqdm


def CreatePolygonz(points):
    """创建面文件：points为一个列表；返回面文件读取信息"""
    outshp = 'temp/temp.shp'
    w = shapefile.Writer(outshp)  # 注意，这里的参数不可以是shapeType = 5，必须是文件路径，否则会报错
    # 设置字段，最大长度为254，C为字符串
    w.field('FIRST_FLD')
    w.field('SECOND_FLD', 'C', '40')
    w.poly([points])
    w.record('First', 'Point')
    # 保存
    w.close()
    # 设置投影，通过.prj文件设置，需要写入一个wkt字符串
    ##gdal的GetProjection()返回的是wkt字符串，需要ImportFromWkt
    # projstr="""PROJCS["WGS_1984_UTM_zone_50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]'"""
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(4326)
    # 或 proj.ImportFromProj4(proj4str)等其他的来源
    wkt = proj.ExportToWkt()
    # 写出prj文件
    f = open(outshp.replace(".shp", ".prj"), 'w')
    f.write(wkt)
    f.close()
    shp = shapefile.Reader(outshp)
    return shp


if __name__ == '__main__':
    # 取得深圳的shape
    sz_shp = shapefile.Reader(r'shp/水体')
    shapes = sz_shp.shapes()
    print('共计{}个面要素'.format(len(shapes)))
    in_shape_points = []
    for index in tqdm(range(len(shapes))):
        geo = shapes[index]
        new_shp = CreatePolygonz(geo.points)
        lon_0, lat_0, lon_1, lat_1 = new_shp.bbox

        linear_lon = np.linspace(lon_0, lon_1, 50)  # 经向50个点
        linear_lat = np.linspace(lat_0, lat_1, 25)  # 纬向25个点
        grid_lon, grid_lat = np.meshgrid(linear_lon, linear_lat)  # 构成了一个坐标矩阵，实际上也就是一个网格，两者是同型矩阵
        flat_lon = grid_lon.flatten()  # 将坐标展成一维
        flat_lat = grid_lat.flatten()

        flat_points = np.column_stack((flat_lon, flat_lat))

        in_shape_points = []
        for pt in flat_points:
            # make a point and see if it's in the polygon
            if geometry.Point(pt).within(geometry.shape(shapes)):
                in_shape_points.append(pt)
        selected_lon = [elem[0] for elem in in_shape_points]
        selected_lat = [elem[1] for elem in in_shape_points]