不想写了
# 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