-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgeo_lab.py
210 lines (182 loc) · 7.95 KB
/
geo_lab.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
from Clawer_Base.logger import logger
from math import radians, cos, sin, asin, sqrt, hypot
import numpy as np
import shapefile
import pandas as pd
class Geo_Point:
def __init__(self, lng, lat):
self.lng = lng
self.lat = lat
def calc_distance(self,another_point):
# 将十进制度转化为弧度
lng1, lat1, lng2, lat2 = map(radians, [self.lng, self.lat, another_point.lng, another_point.lat])
dlng = lng2 - lng1
dlat = lat2 - lat1
h = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlng/2)**2
earth_radius = 6371000
distance = 2 * asin(sqrt(h)) * earth_radius
return distance
def convert_to_rect(self, distance):
# distance 单位为m
lng_per_meter = 0.00001141
lat_per_meter = 0.00000899
lng1 = self.lng - (lng_per_meter * distance)
lat1 = self.lat - (lat_per_meter * distance)
lng2 = self.lng + (lng_per_meter * distance)
lat2 = self.lat + (lat_per_meter * distance)
return Rectangle(lng1, lat1, lng2, lat2)
def convert_to_dict(self, key_name):
a_dict = {}
a_dict[key_name] = str(self.lng)+ ',' + str(self.lat)
return a_dict
class Geo_line:
def __init__(self,lng1, lat1, lng2, lat2):
self.point_1 = Geo_Point(lng1, lat1)
self.point_2 = Geo_Point(lng2, lat2)
self.length = self.point_1.calc_distance(self.point_2)
def convert_to_point(self, accuracy):
divide = int(self.length / accuracy) + 1
lng_min = min([self.point_1.lng, self.point_2.lng])
lng_max = max([self.point_1.lng, self.point_2.lng])
lat_min = min([self.point_1.lat, self.point_2.lat])
lat_max = max([self.point_1.lat, self.point_2.lat])
lng_list = np.linspace(lng_min, lng_max, divide)
lat_list = np.linspace(lat_min, lat_max, divide)
point_list = [Geo_Point(*i) for i in zip(lng_list, lat_list)]
return point_list
class Rectangle:
def __init__(self, lng1=0, lat1=0, lng2=0, lat2=0):
self.lng1 = lng1
self.lat1 = lat1
self.lng2 = lng2
self.lat2 = lat2
self.left_up = Geo_Point(lng1, lat2)
self.left_down = Geo_Point(lng1, lat1)
self.right_up = Geo_Point(lng2, lat2)
self.right_down = Geo_Point(lng2, lat1)
self.center = Geo_Point((lng1+lng2)/2,(lat1+lat2)/2)
self.wide = self.calc_wide()
self.high = self.calc_high()
self.radius = self.out_circle_radius()
def calc_wide(self):
return self.left_up.calc_distance(self.right_up)
def calc_high(self):
return self.left_up.calc_distance(self.left_down)
def out_circle_radius(self):
radius = hypot(self.wide, self.high)/2
return radius
def expand(self, distance, lng1=0, lat1=0, lng2=0, lat2=0):
# distance 单位为m
lng_per_meter = 0.00001141
lat_per_meter = 0.00000899
ex_lng1 = round([(self.lng1 - (lng_per_meter * distance)), lng1][lng1 != 0],6)
ex_lng2 = round([(self.lng2 + (lng_per_meter * distance)), lng2][lng2 != 0],6)
ex_lat1 = round([(self.lat1 - (lat_per_meter * distance)), lat1][lat1 != 0],6)
ex_lat2 = round([(self.lat2 + (lat_per_meter * distance)), lat2][lat2 != 0],6)
logger.info('拓展矩形为(%s,%s,%s,%s)'%(ex_lng1, ex_lat1, ex_lng2, ex_lat2))
return Rectangle(ex_lng1, ex_lat1, ex_lng2, ex_lat2)
def convert_to_lines(self):
self.left_line = Geo_line(self.lng1, self.lat1, self.lng1, self.lat2)
self.right_line = Geo_line(self.lng2, self.lat1, self.lng2, self.lat2)
self.up_line = Geo_line(self.lng1, self.lat2, self.lng2, self.lat2)
self.down_line = Geo_line(self.lng1, self.lat1, self.lng2, self.lat1)
return {'left': self.left_line,
'right': self.right_line,
'up': self.up_line,
'down': self.down_line
}
def convert_to_eight_points(self):
self.middle_up = Geo_Point(self.center.lng,self.lat2)
self.middle_down = Geo_Point(self.center.lng,self.lat1)
self.middle_left = Geo_Point(self.lng1,self.center.lat)
self.middle_right = Geo_Point(self.lng2, self.center.lat)
return [self.left_up,
self.left_down,
self.right_up,
self.right_down,
self.middle_up,
self.middle_down,
self.middle_left,
self.middle_right
]
def convert_to_outline_square(self):
deg_per_meter = (self.lng2 - self.lng1)/self.wide
new_lng1 = self.center.lng - (0.5 * deg_per_meter * self.high)
new_lng2 = self.center.lng + (0.5 * deg_per_meter * self.high)
return Rectangle(new_lng1, self.lat1, new_lng2, self.lat2)
def read_from_shp(self, file_path):
sf = shapefile.Reader(file_path)
self.bbox = sf.bbox
return Rectangle(*tuple(self.bbox))
def divided_into_four(self):
rect_left_down = Rectangle(self.left_down.lng,
self.left_down.lat,
self.center.lng,
self.center.lat)
rect_right_up = Rectangle(self.center.lng,
self.center.lat,
self.right_up.lng,
self.right_up.lat)
rect_left_up = Rectangle(self.left_down.lng,
self.center.lat,
self.center.lng,
self.right_up.lat)
rect_right_down = Rectangle(self.center.lng,
self.left_down.lat,
self.right_up.lng,
self.center.lat)
return [rect_left_down, rect_right_up, rect_left_up, rect_right_down]
def convert_to_param_dict(self):
a_dict = {}
a_dict['location'] = str(self.center.lat) + ',' + str(self.center.lng)
a_dict['radius'] = str(self.radius)
return a_dict
def convert_to_df_dict(self):
a_dict = {}
a_dict['lng'] = self.center.lng
a_dict['lat'] = self.center.lat
a_dict['radius'] = self.radius
return a_dict
class Sample_Generator:
def __init__(self, region_name, category=''):
self.count_sati_rects = []
self.radius_sati_rects = []
self.category = category
self.region_name = region_name
def filter_radius(self, rects, radius):
print('开始生成 %s 采样区域' % self.region_name)
loop_count = 0
while rects:
loop_count += 1
print('第 %s 计算' % loop_count)
rect = rects.pop()
if rect.radius > radius:
rects.extend(rect.divided_into_four())
else:
self.radius_sati_rects.append(rect)
print(u"%s 生成少于 %s 采样点 %s 个" % (self.region_name, radius, len(self.radius_sati_rects)))
def filter_count(self, rects, res_num):
"""根据具体爬虫进行结果筛选"""
print('跳过')
pass
if __name__ == "__main__":
# a_rect = Rectangle().read_from_shp(r'D:\program_lib\GDPOI\guangzhou\广州shapefile')
# square = a_rect.convert_to_outline_square()
# print(square.lng1)
# print(square.lng2)
# print(square.lat1)
# print(square.lat2)
df = pd.read_excel(r'D:\program_lib\BD_scenery\result\青烟威合并.xlsx_group.xlsx')
dict_list = df.to_dict('records')
res_list = []
for i in dict_list:
point1_lng = i['gcj_lng']
point1_lat = i['gcj_lat']
point2_lng = i['gcj_lng1']
point2_lat = i['gcj_lat1']
geo_point = Geo_Point(point1_lng, point1_lat)
geo_point1 = Geo_Point(point2_lng, point2_lat)
i['distance'] = geo_point.calc_distance(geo_point1)
res_list.append(i)
res_df = pd.DataFrame(res_list)
res_df.to_excel(r'D:\program_lib\BD_scenery\result\青烟威合并distance.xlsx')