Huge number of geometric object takes time to load (Especially, loading standardized regional meshes are time consuming). Therefore, I implemented standardized regional mesh generator which generates meshes according to input boundary box.
#!/usr/bin/env python
#coding: utf-8
import sys
import datetime
from rtree import index
from rtree.index import Rtree
from shapely.geometry import Polygon
class SrmGenerator():
def get_mesh1st_to_polygon(self, bbox):
return dict([(mesh1st, Polygon(self._get_rect(*bbox)))
for mesh1st, bbox in self.generate_mesh1st_to_bbox(bbox).items()])
def generate_mesh1st_to_bbox(self, bbox):
d_mesh1st_to_bbox = self._get_d_mesh1st_to_bbox()
rindex = self._create_rtree_index(d_mesh1st_to_bbox)
return dict([(mesh1st, d_mesh1st_to_bbox[mesh1st])
for mesh1st in list(rindex.intersection(bbox))])
def get_mesh2nd_to_polygon(self, bbox):
return dict([(mesh2nd, Polygon(self._get_rect(*bbox)))
for mesh2nd, bbox in self.generate_mesh2nd_to_bbox(bbox).items()])
def generate_mesh2nd_to_bbox(self, bbox):
d_mesh1st_to_bbox = self.generate_mesh1st_to_bbox(bbox)
d_mesh2nd_to_bbox = self._get_d_mesh2nd_to_bbox(d_mesh1st_to_bbox)
rindex = self._create_rtree_index(d_mesh2nd_to_bbox)
return dict([(mesh2nd, d_mesh2nd_to_bbox[mesh2nd])
for mesh2nd in list(rindex.intersection(bbox))])
def get_mesh3rd_to_polygon(self, bbox):
return dict([(mesh3rd, Polygon(self._get_rect(*bbox)))
for mesh3rd, bbox in self.generate_mesh3rd_to_bbox(bbox).items()])
def generate_mesh3rd_to_bbox(self, bbox):
d_mesh2nd_to_bbox = self.generate_mesh2nd_to_bbox(bbox)
d_mesh3rd_to_bbox = self._get_d_mesh3rd_to_bbox(d_mesh2nd_to_bbox)
rindex = self._create_rtree_index(d_mesh3rd_to_bbox)
return dict([(mesh3rd, d_mesh3rd_to_bbox[mesh3rd])
for mesh3rd in list(rindex.intersection(bbox))])
def get_mesh4th_to_polygon(self, bbox):
return dict([(mesh4th, Polygon(self._get_rect(*bbox)))
for mesh4th, bbox in self.generate_mesh4th_to_bbox(bbox).items()])
def generate_mesh4th_to_bbox(self, bbox):
d_mesh3rd_to_bbox = self.generate_mesh3rd_to_bbox(bbox)
d_mesh4th_to_bbox = self._get_d_mesh4th_to_bbox(d_mesh3rd_to_bbox)
rindex = self._create_rtree_index(d_mesh4th_to_bbox)
return dict([(mesh4th, d_mesh4th_to_bbox[mesh4th])
for mesh4th in list(rindex.intersection(bbox))])
def _create_rtree_index(self, d_meshid_to_bbox):
p = index.Property()
rindex = index.Index(properties=p)
[rindex.add(meshid, bbox) for meshid, bbox in d_meshid_to_bbox.items()]
return rindex
def _get_d_mesh1st_to_bbox(self):
l_x = [str(x).zfill(2) for x in range(20, 50)]
l_y = [str(y).zfill(2) for y in range(30, 70)]
return dict([self._generate_mesh1st_bbox_pair(x, y)
for x, y in self._get_combinations(l_x, l_y)])
def _generate_mesh1st_bbox_pair(self, x, y):
id_mesh1st = int("{0}{1}".format(y, x))
left = float("1{0}".format(x))
bottom = float(y) / 1.5
right = float("1{0}".format(str(int(x) + 1)))
top = float(int(y) + 1) / 1.5
return (id_mesh1st, (left, bottom, right, top))
def _get_d_mesh2nd_to_bbox(self, d_mesh1st_to_bbox, num_split=8):
fn = self._get_meshid_for_2nd_3rd
return self._get_d_meshid_to_bbox(d_mesh1st_to_bbox, num_split, fn)
def _get_d_mesh3rd_to_bbox(self, d_mesh2nd_to_bbox, num_split=10):
fn = self._get_meshid_for_2nd_3rd
return self._get_d_meshid_to_bbox(d_mesh2nd_to_bbox, num_split, fn)
def _get_d_mesh4th_to_bbox(self, d_mesh3rd_to_bbox, num_split=2):
fn = self._get_meshid_for_4th
return self._get_d_meshid_to_bbox(d_mesh3rd_to_bbox, num_split, fn)
def _get_d_meshid_to_bbox(self, d_meshcode_to_bbox, num_split, fn):
D = {}
for meshcode, bbox in d_meshcode_to_bbox.items():
left, bottom, right, top = bbox
x_delta = (right - left) / num_split
y_delta = (top - bottom) / num_split
x_coord = [left] + [left+x_delta*i for i in range(1, num_split)] + [right]
y_coord = [bottom] + [bottom+y_delta*i for i in range(1, num_split)] + [top]
for i, j in self._get_combinations(range(len(y_coord)-1), range(len(x_coord)-1)):
new_meshcode = fn(meshcode, i, j)
bbox = (x_coord[j], y_coord[i], x_coord[j+1], y_coord[i+1])
D.update({new_meshcode: bbox})
return D
def _get_meshid_for_2nd_3rd(self, meshcode, i, j):
return int(str(meshcode) + str(i) + str(j))
def _get_meshid_for_4th(self, meshcode, i, j):
return int(str(meshcode) + self._get_mesh4th_last_code(i, j))
def _get_mesh4th_last_code(self, y, x):
if (y == 0 and x == 0):
return '1'
elif (y == 0 and x == 1):
return '2'
elif (y == 1 and x == 0):
return '3'
elif (y == 1 and x == 1):
return '4'
else:
print "unexpected error in _get_mesh4th_last_code"
sys.exit(0)
def _get_combinations(self, L1, L2):
ret = []
for e1 in L1 :
for e2 in L2:
ret.append([e1, e2])
return ret
def _get_rect(self, left, bottom, right, top):
return ((left, bottom),
(left, top),
(right, top),
(right, bottom),
(left, bottom))
if __name__ == "__main__":
L = 130.0
R = 140.0
B = 34.0
T = 34.2
bbox = (L, B, R, T)
from datetime import datetime
s = datetime.now()
print len(SrmGenerator().get_mesh4th_to_polygon(bbox))
e = datetime.now()
print(e-s)
$ python generate_area_srg.py
80100
0:00:07.602852
