python实现创建范围内点的100m正方形格网--地理坐标4326

不想写了

# encoding:  utf8

import math
import arcpy
from arcpy import env


# 创建小方格
# coordinates: 小方格的点集坐标值
# outfile: 小方格图层的名称
def createpolygon(coordinates, outfile):
    # 空间参考
    outSR = arcpy.SpatialReference(3857)

    features = []
    for feature in coordinates:
        # Create a Polygon object based on the array of points
        # Append to the list of Polygon objects
        features.append(
            arcpy.Polygon(
                arcpy.Array([arcpy.Point(*coords) for coords in feature])))
    # 保存要素到工作空间
    arcpy.CopyFeatures_management(features, outfile)
    # 定义投影
    arcpy.DefineProjection_management(outfile, outSR)


# 4326 to web mercator

def lonLatToMercator(lon, lat):
    mercator = {'x': 0, 'y': 0}
    earthRad = 6378137.0
    mercator['x'] = lon * math.pi / 180 * earthRad
    a = lat * math.pi / 180
    mercator['y'] = earthRad / 2 * math.log((1.0 + math.sin(a)) / (1.0 - math.sin(a)))
    return mercator


# --------------------------------

# 工作空间

environment = u'D:/学习/提取小方格/data4326.gdb'
env.workspace = environment

outfile = "envelope4326"

dataset = 'range4326'
extent = arcpy.Describe(dataset).extent
# 获取范围要素的坐标值
# Xmin, Ymin, Xmax, Ymax
xy = lonLatToMercator(extent.XMax, extent.YMax)
Xmax = xy['x']
Ymax = xy['y']

# 读取点图层中的点要素获取坐标
pointfc = 'point4326'
coors = []
with arcpy.da.SearchCursor(pointfc, 'SHAPE@XY') as cursor:
    for row in cursor:
        print row
        tempgeo = lonLatToMercator(row[0][0], row[0][1])
        Xp = tempgeo['x']
        Yp = tempgeo['y']

        LeftTop = [((Xmax - Xp) % 100) + Xp - 100, ((Ymax - Yp) % 100) + Yp]
        RightTop = [((Xmax - Xp) % 100) + Xp, ((Ymax - Yp) % 100) + Yp]
        RightBottom = [((Xmax - Xp) % 100) + Xp - 100, ((Ymax - Yp) % 100) + Yp - 100]
        LeftBottom = [((Xmax - Xp) % 100) + Xp, ((Ymax - Yp) % 100) + Yp - 100]

        tempcoor = [LeftTop, RightTop, LeftBottom, RightBottom]
        coors.append(tempcoor)

print coors

# coor = [[[1.1, 2.1], [2.1, 4.1], [3.1, 7.1]],
#                 [[6.2, 8.2], [5.2, 7.2], [7.2, 2.2], [9.2, 5.2]]]
#
# coor = [[[11000604.857985009, 3141841.7007190087], [11000704.857985009, 3141841.7007190087],
#         [11000704.857985009, 3141741.7007190087], [11000604.857985009, 3141741.7007190087]]]
createpolygon(coors, outfile)



# 这个没用上，3857的坐标，这里需要改椭球体参数，就用了上面那个
'''
2000国家大地坐标系的大地测量基本常数分别为:
长半轴 a = 6378137 m        
地球引力常数 GM =3.986004418×1014m3s-2
扁率f = 1/ 298.257222101 = 0.0033528106811823  
地球自转角速度X =7.292115×10-5rad s-1
'''


def latlng2xy(latitude, longitude):
    a = 6378137.0
    e2 = 0.00669342162296594323
    # 将经纬度转换为弧度
    latitude2Rad = (math.pi / 180.0) * latitude

    beltNo = int((longitude + 1.5) / 3.0)  # 计算3度带投影度带号
    L = beltNo * 3  # 计算中央经线
    l0 = longitude - L  # 经差
    tsin = math.sin(latitude2Rad)
    tcos = math.cos(latitude2Rad)
    t = math.tan(latitude2Rad)
    m = (math.pi / 180.0) * l0 * tcos
    et2 = e2 * pow(tcos, 2)
    et3 = e2 * pow(tsin, 2)
    X = 111132.9558 * latitude - 16038.6496 * math.sin(2 * latitude2Rad) + 16.8607 * math.sin(
        4 * latitude2Rad) - 0.0220 * math.sin(6 * latitude2Rad)
    N = a / math.sqrt(1 - et3)

    x = X + N * t * (0.5 * pow(m, 2) + (5.0 - pow(t, 2) + 9.0 * et2 + 4 * pow(et2, 2)) * pow(m, 4) / 24.0 + (
            61.0 - 58.0 * pow(t, 2) + pow(t, 4)) * pow(m, 6) / 720.0)
    y = 500000 + N * (m + (1.0 - pow(t, 2) + et2) * pow(m, 3) / 6.0 + (
            5.0 - 18.0 * pow(t, 2) + pow(t, 4) + 14.0 * et2 - 58.0 * et2 * pow(t, 2)) * pow(m, 5) / 120.0)

    return x, y

转坐标有点烦，最开始用常数去乘误差太大，要不得，后面问了大佬，把下面这个改成py的就行了，原来是js的。

参考：https://blog.csdn.net/semian7633/article/details/109258717