moldovanGrid.py

import shapely.geometry
import shapely.affinity
import geopandas as gpd

origin = shapely.geometry.Point(67596.000000, 36694.000000)
pixel_count = 256
units_per_pixel_for_each_zoom_level = [
    352.77758727788068426889,
    176.38879363894034213445,
    88.19439681947017106722,
    35.27775872778806132146,
    17.63887936389403066073,
    8.81943968194701533037,
    3.52777587277880622096,
    1.76388793638940311048,
    0.70555517455576133301,
    0.35277758727788066651,
    0.17638879363894033325,
    0.08819439699999999382,
    0.03527775880000000308,
]  # from https://geoportal.md/tiles/Tipulsolului/

for z, units_per_pixel in enumerate(units_per_pixel_for_each_zoom_level[:4]):
    ps = []
    ids = []
    side_length = units_per_pixel * pixel_count

    p = shapely.geometry.box(
        origin.x,
        origin.y,
        origin.x + (units_per_pixel * pixel_count),
        origin.y + (units_per_pixel * pixel_count),
    )

    for i in range(0, 8*(z+1)):
        for j in range(0, 8*(z+1)):
            ps.append(
                shapely.affinity.translate(p, xoff=i * side_length, yoff=j * side_length)
            )
            ids.append(f"{z}/{i}/{j}")

    gpd.GeoDataFrame.from_dict({"geometry": ps, "id": ids}).set_crs(4026).to_file(
        f"grid_z{z}.geojson"
    )