forked from gsi-cyberjapan/vector_tiles_convert
-
Notifications
You must be signed in to change notification settings - Fork 0
/
clipping.py
204 lines (180 loc) · 6.4 KB
/
clipping.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
from os.path import join,exists,basename,splitext
import math
import pyproj
import json
import geojson
import globalmaptiles
from load import loader
from shapely.geometry import mapping,box
from datetime import datetime as dt
from utils import u_makedirs as makedirs
def conv_latlon(lat,lon,inv=None):
if inv is None:
src_proj = pyproj.Proj(init="epsg:4612")
dst_proj = pyproj.Proj(init="epsg:3857")
else:
dst_proj = pyproj.Proj(init="epsg:4612")
src_proj = pyproj.Proj(init="epsg:3857")
nlon,nlat = pyproj.transform(src_proj,dst_proj,lon,lat)
return (nlat,nlon)
# get latlon on bottom-left corner of a mesh
def get_latlon(code):
if len(code) < 4:
print("code is too short!")
lat = float(code[0:2]) / 1.5
lon = float(code[2:4]) + 100.0
if len(code) >= 6:
lat += float(code[4]) * 5.0 / 60.0
lon += float(code[5]) * 7.5 / 60.0
if len(code) >= 8:
lat += float(code[6]) * 30.0 / 3600.0
lon += float(code[7]) * 45.0 / 3600.0
if len(code) >= 9:
if code[8] == "2" or code[8] == "4":
lon += 22.5 / 3600.0
if code[8] == "3" or code[8] == "4":
lat += 15.0 / 3600.0
return (lat,lon)
#get latlon on all corners
def get_latlon_corners(code):
minlat,minlon = get_latlon(code)
if len(code) == 4:
maxlat = minlat + 40.0/60.0
maxlon = minlon + 1.0
elif len(code) == 6:
maxlat = minlat + 5.0/60.0
maxlon = minlon + 7.5/60.0
elif len(code) == 8:
maxlat = minlat + 30.0/3600.0
maxlon = minlon + 45.0/3600.0
elif len(code) == 9:
maxlat = minlat + 15.0/3600.0
maxlon = minlon + 22.5/3600.0
return (minlat,maxlat,minlon,maxlon)
def get_tile_id_corner(code,zoom):
if not isinstance(zoom,int):
zoom=int(zoom)
minlat,maxlat,minlon,maxlon = get_latlon_corners(code)
min_my,min_mx = conv_latlon(minlat,minlon)
max_my,max_mx = conv_latlon(maxlat,maxlon)
mercator = globalmaptiles.GlobalMercator()
min_tx,min_ty = mercator.MetersToTile(min_mx,min_my,zoom)
max_tx,max_ty = mercator.MetersToTile(max_mx,max_my,zoom)
return [[min_tx,max_tx],[min_ty,max_ty]]
def get_clip_bound(x,y,z):
mercator = globalmaptiles.GlobalMercator()
bounds = mercator.TileBounds( x, y, z)
return bounds
def conv_proj(obj,coord_tab):
ret,coord_tab = rec_in_conv_proj(obj['geometry']['coordinates'],coord_tab)
obj['geometry']['coordinates'] = ret
return obj,coord_tab
def in_conv_proj(obj,coord_tab,inv=None):
if not isinstance(obj,list):
if obj in coord_tab:
return coord_tab[obj]
y,x = conv_latlon(obj[1],obj[0],inv=inv)
return (x,y)
def rec_in_conv_proj(obj,coord_tab,inv=None):
if isinstance(obj[0],float):
res = in_conv_proj(obj,coord_tab,inv)
if inv is None and not res in coord_tab:
coord_tab[res] = (obj[0],obj[1])
return res,coord_tab
elif not obj is None:
ret = []
for o in obj:
z,coord_tab = rec_in_conv_proj(o,coord_tab,inv)
ret.append(z)
return ret,coord_tab
def inv_conv_proj(obj,coord_tab):
newobj = {'geometry':{}}
newobj['properties'] = obj['properties']
newobj['type'] = obj['type']
newobj['geometry']['type'] = obj['geometry']['type']
newobj['geometry']['coordinates'],coord_tab = rec_in_conv_proj(obj['geometry']['coordinates'],coord_tab,1)
return newobj
def json_dump(dic,out,dpath,z,x,y):
mercator = globalmaptiles.GlobalMercator()
idx,idy = mercator.GoogleTile(x,y,z)
l = None
for k,v in dic.items():
if not v:
continue
path = join(out,dpath,k,str(z),str(idx))
makedirs(path)
json_dic = {'type':'FeatureCollection', 'features':v}
txt = json.dumps(json_dic,indent=2,ensure_ascii=False)
p = join(path,str(idy)+'.geojson')
fw = open(p,'w')
fw.write(txt.encode('utf_8'))
fw.close()
def main_shp(zoom,path,out,features=None,test=False):
assert( path.endswiths(".shp"))
geopath = path[:-3] + "geojson"
cmd = "ogr2ogr -f GeoJSON " + geopath + " " + path
system(cmd)
gobj = geojson.load(geopath)
def main(zoom,path,out,features=None,test=False,code=None,shpfile=False):
start = dt.now()
a = loader(path,shpfile)
print "TIME LOAD:",dt.now() - start
if code is None:
code = str(a.code)
trng = get_tile_id_corner(code,zoom)
# check featurelist
for x in features:
if x not in a.gobj.keys():
print("X",x)
print("FS",sorted(features))
print("KS",sorted(a.gobj.keys()))
assert(x in a.gobj.keys())
# convert coordinates
coord_tab = {}
for f in a.gobj:
if not features or f in features:
for i in range(len(a.gobj[f])):
a.gobj[f][i],coord_tab = conv_proj(a.gobj[f][i],coord_tab)
# generate index
idx = {}
vals = {}
for f in a.gobj:
if not features or f in features:
idx[f] = a.rindex(f,rtree=not test)
vals[f] = a.gobj[f]
print "TIME INDEX:",dt.now() - start
# generate bounds
bs = []
bxs = []
for xid in range(trng[0][0],trng[0][1]+1):
for yid in range(trng[1][0],trng[1][1]+1):
bounds = get_clip_bound(xid,yid,int(zoom))
bs.append((xid,yid,tuple(bounds)))
bxs.append((xid,yid,box(*bounds)))
res = {}
if test:
for f in vals:
#print "EXTRACT",f,
tmp = a.extract3(bxs,vals[f])
#print "DONE"
for xid,yid in tmp.keys():
res[xid,yid,f] = tmp[xid,yid]
else:
for xid,yid,bounds in bs:
res[xid,yid] = {}
for f in vals:
#print "EXTRACT",bounds,f,
res[xid,yid][f] = a.extract2(bounds,vals[f],idx[f])
#print "DONE"
print "TIME EXTRACT:",dt.now() - start
for k in res:
for kk in res[k]:
res[k][kk] = [inv_conv_proj(x,coord_tab) for x in res[k][kk]]
print "TIME COORD:",dt.now() - start
for xid in range(trng[0][0],trng[0][1]+1):
for yid in range(trng[1][0],trng[1][1]+1):
dumpres = res[xid,yid]
if dumpres:
json_dump(dumpres,out,code,int(zoom),xid,yid)
print "TIME DUMP:",dt.now() - start
return code