GISTOOLS.md

GISTOOLS

install

python3 -m pip uninstall gistools -y 
python3 -m pip install git+https://gist.github.com/tomowarkar/598768818ed7d2e6ca94f4011dc065a6

example

from gistools import *


shinjuku = (35.70078, 139.71475)
tsukuba = (36.10377, 140.08786)

shinjuku_mesh = coord2mesh(*shinjuku)
tsukuba_coord, _ = mesh2coord(coord2mesh(*tsukuba))
dist = haversine_distance(shinjuku, tsukuba)
theta = math.degrees(azimuth(shinjuku, tsukuba))

print(
    "新宿 基準地域メッシュ:\t%s\n筑波 緯度経度:\t%s\n距離 [km]:\t%.2f\n角度 [°]:\t%.2f\n"
    % (shinjuku_mesh, tsukuba_coord, dist, theta)
)

新宿 基準地域メッシュ:  53394547
筑波 緯度経度:  (36.1, 140.0875)
距離 [km]:      56.07
角度 [°]:       36.76

def calc_real_scale(
    code: str, distance
):
    mesh = Mesh.from_code(code)
    lat, lng = mesh.latitude, mesh.longitude
    h, w = mesh.height, mesh.width

    width = distance((lat, lng), (lat, lng + w))
    height = distance((lat, lng), (lat + h, lng))
    return width, height

w, h = calc_real_scale("39272554", lambert_andoyer_distance)
print("那覇\n横 [km]:\t%.3f\n縦 [km]:\t%.3f\n面積 [km^2]:\t%.3f" % (w, h, w * h))

w, h = calc_real_scale("64414277", lambert_andoyer_distance)
print("札幌\n横 [km]:\t%.3f\n縦 [km]:\t%.3f\n面積 [km^2]:\t%.3f" % (w, h, w * h))

那覇
横 [km]:	1.249
縦 [km]:	0.923
面積 [km^2]:	1.153
札幌
横 [km]:	1.018
縦 [km]:	0.926
面積 [km^2]:	0.943

gistools.py

from __future__ import annotations

import math
from dataclasses import dataclass
from typing import Callable


@dataclass
class Coordinate:
    latitude: float
    longitude: float

    def __iter__(self):
        return iter((self.latitude, self.longitude))


@dataclass
class Mesh(Coordinate):
    code: str
    scale: int = 1

    def __iter__(self):
        return iter(self.center)

    @classmethod
    def from_code(cls, code: str):
        (lat, lng), scale = mesh2coord(code)
        return cls(latitude=lat, longitude=lng, code=code, scale=scale)

    @property
    def width(self):
        return 1 / self.scale

    @property
    def height(self):
        return 2 / 3 / self.scale

    @property
    def center(self):
        return self.latitude + self.height / 2, self.longitude + self.width / 2

    def coordinates(self, lnglat: bool = True):
        s = slice(None, None, 1 if lnglat else -1)
        return [
            [self.longitude, self.latitude][s],
            [self.longitude + self.width, self.latitude][s],
            [self.longitude + self.width, self.latitude + self.height][s],
            [self.longitude, self.latitude + self.height][s],
            [self.longitude, self.latitude][s],
        ]


def coord2mesh(latitude: float, longitude: float, key: int = 3):
    a, _ = divmod(latitude * 60, 40)
    c, _ = divmod(_, 5)
    e, _ = divmod(_ * 60, 30)
    i, _ = divmod(_, 30 / 2)
    k, _ = divmod(_, 30 / 4)
    m, _ = divmod(_, 30 / 8)

    b, _ = divmod(longitude - 100, 1)
    d, _ = divmod(_ * 60, 7.5)
    f, _ = divmod(_ * 60, 45)
    j, _ = divmod(_, 45 / 2)
    l, _ = divmod(_, 45 / 4)
    n, _ = divmod(_, 45 / 8)

    # fmt:off
    code = "%d%d%d%d%d%d%d%d%d" % (a, b, c, d, e, f, (i * 2) + (j + 1), (k * 2) + (l + 1), (m * 2) + (n + 1))
    # fmt: on
    thd = {
        1: 4,  # 80km
        2: 6,  # 10km
        3: 8,  # 　1km
        4: 9,  # 　500m
        5: 10,  # 250m
        6: 11,  # 125m
    }.get(key, None)
    return code[:thd]


# def coord2mesh(latitude: float, longitude: float, key: int = 3):
#     lat, lng = latitude, longitude
#     code = ""

#     if 0 < key:  # 80km　メッシュ
#         a, lat = divmod(lat * 60, 40)
#         b, lng = divmod(lng - 100, 1)
#         code += "%d%d" % (a, b)
#     if 1 < key:  # 10km　メッシュ
#         a, lat = divmod(lat, 5)
#         b, lng = divmod(lng * 60, 7.5)
#         code += "%d%d" % (a, b)
#     if 2 < key:  # 1km　メッシュ
#         a, lat = divmod(lat * 60, 30)
#         b, lng = divmod(lng * 60, 45)
#         code += "%d%d" % (a, b)

#     for i in range(key - 3):  # 1/2**(1+i) km　メッシュ
#         scale = 2 ** (1 + i)
#         a, lat = divmod(lat, 30 / scale)
#         b, lng = divmod(lng, 45 / scale)
#         code += "%d" % ((a * 2) + (b + 1))
#     return code, (lat, lng)


def mesh2coord(code: str, center: bool = False):
    lat, lng, scale = 0, 100, 1
    len_ = len(code)

    if len_ < 4:
        raise ValueError
    if 3 < len_:
        lat += int(code[0:2]) * 2 / 3 / scale
        lng += int(code[2:4]) / scale
    if 5 < len_:
        scale *= 8
        lat += int(code[4:5]) * 2 / 3 / scale
        lng += int(code[5:6]) / scale
    if 7 < len_:
        scale *= 10
        lat += int(code[6:7]) * 2 / 3 / scale
        lng += int(code[7:8]) / scale
    for rune in code[8:]:
        scale *= 2
        lat += (int(rune) in (3, 4)) * 2 / 3 / scale
        lng += (int(rune) in (2, 4)) / scale
    if center:
        lat += 1 / 3 / scale
        lng += 1 / scale / 2

    return (lat, lng), scale


def azimuth(from_: tuple[float, float], to_: tuple[float, float]):
    """
    from_ : (lat, lng)
    to_ : (lat, lng)

    returns : radian
    """
    (y1, x1), (y2, x2) = map(math.radians, from_), map(math.radians, to_)
    dx = x2 - x1
    s, c = math.sin, math.cos
    theta = math.atan2(c(y1) * math.tan(y2) - s(y1) * c(dx), s(dx))
    return math.pi / 2 - theta


def haversine_distance(from_: tuple[float, float], to_: tuple[float, float]):
    """
    from_ : (lat, lng)
    to_ : (lat, lng)

    returns : distance [km]
    """
    (y1, x1), (y2, x2) = map(math.radians, from_), map(math.radians, to_)
    dx, dy = x2 - x1, y2 - y1
    s, c = math.sin, math.cos

    dist = math.asin(math.sqrt(s(dy / 2) ** 2 + c(y1) * c(y2) * s(dx / 2) ** 2))
    r = 6378.137
    return dist * 2 * r


def lambert_andoyer_distance(from_: tuple[float, float], to_: tuple[float, float]):
    """
    from_ : (lat, lng)
    to_ : (lat, lng)

    returns : distance [km]
    """
    (y1, x1), (y2, x2) = map(math.radians, from_), map(math.radians, to_)

    a = 6378.140
    b = 6356.755
    f = 1 - b / a
    s, c = math.sin, math.cos

    pa = math.atan(b / a * math.tan(y1))
    pb = math.atan(b / a * math.tan(y2))
    xx = math.acos(s(pa) * s(pb) + c(pa) * c(pb) * c(x1 - x2))
    c1 = (s(xx) - xx) * (s(pa) + s(pb)) ** 2 / c(xx / 2) ** 2
    c2 = (s(xx) + xx) * (s(pa) - s(pb)) ** 2 / s(xx / 2) ** 2
    dr = f / 8 * (c1 - c2)
    return a * (xx + dr)


if __name__ == "__main__":
    shinjuku = (35.70078, 139.71475)
    tsukuba = (36.10377, 140.08786)

    shinjuku_mesh = coord2mesh(*shinjuku)
    assert shinjuku_mesh == "53394547"
    tsukuba_coord, _ = mesh2coord(coord2mesh(*tsukuba))
    assert tsukuba_coord == (36.1, 140.0875)
    dist = haversine_distance(shinjuku, tsukuba)
    assert round(dist, 2) == 56.07
    theta = math.degrees(azimuth(shinjuku, tsukuba))
    assert round(theta, 2) == 36.76
    assert round(math.degrees(azimuth(tsukuba, shinjuku)), 2) == 216.98

    print(
        "新宿 基準地域メッシュ:\t%s\n筑波 緯度経度:\t%s\n距離 [km]:\t%.2f\n角度 [°]:\t%.2f\n"
        % (shinjuku_mesh, tsukuba_coord, dist, theta)
    )
    # 新宿 基準地域メッシュ:  53394547
    # 筑波 緯度経度:  (36.1, 140.0875)
    # 距離 [km]:      56.07
    # 角度 [°]:       36.76

    def calc_real_scale(
        code: str, distance: Callable[[tuple[float, float], tuple[float, float]], float]
    ):
        mesh = Mesh.from_code(code)
        lat, lng = mesh.latitude, mesh.longitude
        h, w = mesh.height, mesh.width

        width = distance((lat, lng), (lat, lng + w))
        height = distance((lat, lng), (lat + h, lng))
        return width, height

    w, h = calc_real_scale("39272554", lambert_andoyer_distance)
    assert (round(w, 3) == 1.249) & (round(h, 3) == 0.923)
    print("那覇\n横 [km]:\t%.3f\n縦 [km]:\t%.3f\n面積 [km^2]:\t%.3f" % (w, h, w * h))

    w, h = calc_real_scale("64414277", lambert_andoyer_distance)
    assert (round(w, 3) == 1.018) & (round(h, 3) == 0.926)
    print("札幌\n横 [km]:\t%.3f\n縦 [km]:\t%.3f\n面積 [km^2]:\t%.3f" % (w, h, w * h))

    # 那覇
    # 横 [km]:	1.249
    # 縦 [km]:	0.923
    # 面積 [km^2]:	1.153
    # 札幌
    # 横 [km]:	1.018
    # 縦 [km]:	0.926
    # 面積 [km^2]:	0.943

gistoolsv2.py

import re
import warnings
from typing import Tuple

Latitude = float
Longitude = float
Coordinate = Tuple[Latitude, Longitude]


def valid_coordinate_jp(p: Coordinate):
    lat, lon = p
    lat_is_valid = 20 <= lat <= 46
    lon_is_valid = 122 <= lon <= 154
    return lat_is_valid & lon_is_valid


def _coord2code(p: Coordinate, key="1km"):
    self = _coord2code
    if key == "80km":
        lat, lon = p
        a, c = divmod(lat * 60, 40)
        b, d = divmod(lon - 100, 1)
        return (a, b), (c, d)
    elif key == "10km":
        codes, (c, d) = self(p, key="80km")
        a, c = divmod(c, 5)
        b, d = divmod(d * 60, 7.5)
        codes += (a, b)
        return codes, (c, d)
    elif key == "1km":
        codes, (c, d) = self(p, key="10km")
        a, c = divmod(c * 60, 30)
        b, d = divmod(d * 60, 45)
        codes += (a, b)
        return codes, (c, d)
    elif key == "500m":
        codes, (c, d) = self(p, key="1km")
        a, c = divmod(c, 15)
        b, d = divmod(d, 22.5)
        codes += ((a * 2) + (b + 1),)
        return codes, (c, d)
    elif key == "250m":
        codes, (c, d) = self(p, key="500m")
        a, c = divmod(c, 7.5)
        b, d = divmod(d, 11.25)
        codes += ((a * 2) + (b + 1),)
        return codes, (c, d)
    elif key == "125m":
        codes, (c, d) = self(p, key="250m")
        a, c = divmod(c, 3.75)
        b, d = divmod(d, 5.625)
        codes += ((a * 2) + (b + 1),)
        return codes, (c, d)
    elif key == "5km":
        codes, (c, d) = self(p, key="10km")
        a, c = divmod(c, 2.5)
        b, d = divmod(d, 3.75)
        codes += ((a * 2) + (b + 1),)
        return codes, (c, d)
    elif key == "2km":
        codes, (c, d) = self(p, key="10km")
        a, c = divmod(c, 1)
        b, d = divmod(d, 1.5)
        codes += (a * 2, b * 2, 5)
        return codes, (c, d)
    else:
        raise ValueError("key: %s is not defined." % key)


def coord2code(p: Coordinate, key="1km"):
    codes, _ = _coord2code(p=p, key=key)
    return "".join(map(str, map(int, codes)))


def to_mesh(latitude: float, longitude: float):
    a, c = divmod(latitude * 60, 40)
    b, d = divmod(longitude - 100, 1)
    e, g = divmod(c, 5)
    f, h = divmod(d * 60, 7.5)
    i, _ = divmod(g * 60, 30)
    j, _ = divmod(h * 60, 45)
    return "".join(map(str, map(int, (a, b, e, f, i, j))))


class MeshCode:
    regex = r"~"
    parent = None
    scale = 1
    
    def __init__(self, code: str):
        self.code = code
        assert self.match(self.code)
        self.latitude, self.longitude = self.code2coord(code)
        
    @property
    def width(self):
        return 1 / self.scale

    @property
    def height(self):
        return 2 / 3 / self.scale

    @property
    def center(self):
        return self.latitude + self.height / 2, self.longitude + self.width / 2

    def coordinates(self):
        return [
            [self.latitude, self.longitude],
            [self.latitude, self.longitude + self.width],
            [self.latitude + self.height, self.longitude + self.width],
            [self.latitude + self.height, self.longitude],
            [self.latitude, self.longitude],
        ]

    @classmethod
    def match(cls, code: str):
        return re.match(cls.regex, code) is not None

    @classmethod
    def _coord2code(cls, codes:Tuple[float], lat:Latitude, lon:Longitude):
        return (), (lat, lon - 100)

    @classmethod
    def _coord2code_(cls, p:Coordinate):
        if cls.parent is None:
            codes, (lat, lon) = cls._coord2code((), *p)
        else:
            codes, (lat, lon) = cls.parent._coord2code_(p)
            codes, (lat, lon) = cls._coord2code(codes, lat, lon)
        return codes, (lat, lon)

    @classmethod
    def coord2code(cls, p:Coordinate):
        codes, _ = cls._coord2code_(p)
        code = "".join(map(str, map(int, codes)))
        return code

    @classmethod
    def from_coord(cls, p:Coordinate):
        code = cls.coord2code(p)
        return cls(code)

    @classmethod
    def _code2coord(cls, code: str):
        return (0, 100)

    @classmethod
    def code2coord(cls, code: str):
        if cls.parent is None:
            return cls._code2coord(code)
        plat, plon = cls.parent.code2coord(code)
        lat, lon = cls._code2coord(code)
        return (lat + plat, lon + plon)


class MeshCode80km(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])$"
    parent = MeshCode
    scale = 1

    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat * 60, 40)
        b, d = divmod(lon, 1)
        codes += (a, b)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code:str):
        lat, lon = map(float, (code[0:2], code[2:4]))
        return lat * 2 / 3, lon


class MeshCode10km(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}$"
    parent = MeshCode80km
    scale = 8
    
    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat, 5)
        b, d = divmod(lon * 60, 7.5)
        codes += (a, b)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = map(float, (code[4:5], code[5:6]))
        return lat * 2 / 3 / 8, lon / 8


class MeshCode5km(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}[1-4]$"
    parent = MeshCode10km
    scale = 16
    
    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat, 2.5)
        b, d = divmod(lon, 3.75)
        codes += ((a * 2) + (b + 1),)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = divmod(int(code[6:7]) - 1, 2)
        return lat * 2 / 3 / 16, lon / 16


class MeshCode2km(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}[02468]{2}5$"
    parent = MeshCode10km
    scale = 40

    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat, 1)
        b, d = divmod(lon, 1.5)
        codes += (a * 2, b * 2, 5)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = map(float, (code[6:7], code[7:8]))
        return lat * 2 / 3 / 80, lon / 80


class MeshCode1km(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}\d{2}$"
    parent = MeshCode10km
    scale = 80

    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat * 60, 30)
        b, d = divmod(lon * 60, 45)
        codes += (a, b)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = map(float, (code[6:7], code[7:8]))
        return lat * 2 / 3 / 80, lon / 80


class MeshCode500m(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}\d{2}[1-4]$"
    parent = MeshCode1km
    scale = 160

    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat, 15)
        b, d = divmod(lon, 22.5)
        codes += ((a * 2) + (b + 1),)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = divmod(int(code[8:9]) - 1, 2)
        return lat * 2 / 3 / 160, lon / 160


class MeshCode250m(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}\d{2}[1-4]{2}$"
    parent = MeshCode500m
    scale = 320

    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat, 7.5)
        b, d = divmod(lon, 11.25)
        codes += ((a * 2) + (b + 1),)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = divmod(int(code[9:10]) - 1, 2)
        return lat * 2 / 3 / 320, lon / 320


class MeshCode125m(MeshCode):
    regex = r"^([3-6][0-9])(2[2-9]|[3-4][0-9]|5[0-5])[0-7]{2}\d{2}[1-4]{3}$"
    parent = MeshCode250m
    scale = 640

    @classmethod
    def _coord2code(cls, codes, lat, lon):
        a, c = divmod(lat, 3.75)
        b, d = divmod(lon, 5.625)
        codes += ((a * 2) + (b + 1),)
        return codes, (round(c, 10), round(d, 10))

    @classmethod
    def _code2coord(cls, code):
        lat, lon = divmod(int(code[10:11]) - 1, 2)
        return lat * 2 / 3 / 640, lon / 640


if __name__ == "__main__":
    import pandas as pd

    # fmt: off
    df = pd.DataFrame(
        [
            ["稚内", 45.41715379, 141.67690228, "6841", "684115", "68411504", "684115041", "6841150411", "68411504112", "6841151", "684115045"],
            ["沖ノ鳥島", 20.425550297436146, 136.081078346884, "3036", "303650", "30365016", "303650161", "3036501612", "30365016122", "3036502", "303650065"],
            ["南鳥島", 24.286699972409636, 153.98078108610335, "3653", "365337", "36533748", "365337481", "3653374814", "36533748144", "3653372", "365337485"],
            ["与那国島", 24.450187920552747, 122.93420216940896, "3622", "362257", "36225744", "362257442", "3622574421", "36225744212", "3622571", "362257445"],
            ["福岡", 33.58972334954021, 130.4207522946408, "5030", "503033", "50303303", "503033034", "5030330343", "50303303432", "5030331", "503033025"],
            ["大阪", 34.702475212431786, 135.4959362881235, "5235", "523503", "52350349", "523503492", "5235034923", "52350349232", "5235032", "523503485"],
            ["東京", 35.68121852142085, 139.76670013272258, "5339", "533946", "53394611", "533946113", "5339461132", "53394611323", "5339461", "533946005"],
            ["青森", 40.828881020977185, 140.73419079481246, "6140", "614015", "61401598", "614015982", "6140159823", "61401598234", "6140154", "614015885"],
            ["札幌", 43.06866073331351, 141.3507799065603, "6441", "644142", "64414288", "644142881", "6441428811", "64414288113", "6441424", "644142885"],
            ["金沢", 36.57804475344767, 136.64795535242703, "5436", "543665", "54366591", "543665912", "5436659124", "54366591241", "5436653", "543665805"],
            ["浮動小数Error", 20.425, 136.075, "3036", "303650", "30365016", "303650161", "3036501611", "30365016111", "3036502", "303650065"],
        ],
        columns=["place", "lat", "lon", "80km", "10km", "1km", "500m", "250m", "125m", "5km", "2km"],
    )

    for _, row in df.iterrows():
        for key, m in zip(
            ["80km", "10km", "1km", "500m", "250m", "125m", "5km", "2km"],
            [MeshCode80km, MeshCode10km, MeshCode1km, MeshCode500m, MeshCode250m, MeshCode125m, MeshCode5km, MeshCode2km]
        ):
        # fmt: on
            p0 = row[['lat', 'lon']]
            code = m.coord2code(p0)
            assert row[key] == code, '%s key="%s": expected="%s" got="%s"' % (row['place'], key, row[key], code)
            
            p1 = m.code2coord(code)
            code2 = m.coord2code(p1)
            p2 = m.code2coord(code2)
            assert code == code2
            assert p1== p2

setup.py

from setuptools import setup


def main():
    setup(
        name="gistools",
        version="1.0.0",
        py_modules=["gistools"],
    )

if __name__ == "__main__":
    main()