geometry and python tips: automatic generation of Japanese standardized regional meshes

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