Shapefileを作るPythonプログラムを無償提供

大量の地物数のShapefileを作成する機会があったので作成したプログラムをご紹介します。
サイズの小さいものや大きなもの地物数等いろいろな条件でShapefileが作成できるので、テスト用にご活用ください。
Shapefileは使用するソフトウェアや環境によって扱えるファイルサイズが異なる場合がありますが、一般的に使われている2GBを超えるShapefileも作成できることを確認しています。

プログラムの仕様

・Shapefileを作成する範囲を緯度経度で指定できる
・PolygonまたはPointのShapefileどちらをつくるか指定できる
・作成する地物数を指定できる
・作成する属性数を指定できる
・Polygonの場合1Polygonの頂点数を指定できる
・文字コードはUTF-8で作成する

プログラムの使い方

プログラムに指定する引数は下記の通りです。
引数１：作成する範囲左上の緯度
引数２：作成する範囲左上の経度
引数３：作成する範囲右下の緯度
引数４：作成する範囲右下の経度
引数５：作成する地物数
引数６：出力ファイルパス
引数７：作成する地物の種類(polygon / point)
オプション：–atr_count：作成する属性数(0～256)
オプション：–num_points：Polygonの頂点数(4以上)
※Lineは作成できません。

コマンド例

いくつかコマンドの例を記載しておきます。

地物数33081、属性情報256項目、地物当たりの頂点数4のPolygonを作る場合
python create_shapefile.py 35.0 135.0 34.0 136.0 33081 ./polygon_33081_256_4.shp polygon --atr_count 256 --num_points 4

地物数84910、属性情報100項目、地物当たりの頂点数4のPolygonを作る場合
python create_shapefile.py 35.0 135.0 34.0 136.0 84910 ./polygon_84910_100_4.shp polygon --atr_count 100 --num_points 4

地物数20、属性情報0項目、地物当たりの頂点数100のPolygonを作る場合
python create_shapefile.py 35.0 135.0 34.0 136.0 20 ./polygon_20_0_100.shp polygon --atr_count 0 --num_points 100

地物数33081、属性情報256項目のPointを作る場合
python create_shapefile.py 35.0 135.0 34.0 136.0 33081 ./point_33081_256.shp point --atr_count 256

地物数84910、属性情報100項目のPointを作る場合
python create_shapefile.py 35.0 135.0 34.0 136.0 84910 ./point_84910_100.shp point --atr_count 100

Pythonプログラム

ファイル名「create_shapefile.py」で保存してご利用ください。

import argparse
import random
import os
from osgeo import ogr, osr
import math
import random
import string

def generate_random_string(length):
    # 使用する文字のセット
    characters = string.ascii_letters + string.digits
    # 指定された長さのランダムな文字列を生成
    random_string = ''.join(random.choice(characters) for _ in range(length))
    return random_string

def generate_random_string(length):
    characters = string.ascii_letters + string.digits
    return ''.join(random.choice(characters) for _ in range(length))

def create_shapefile(left_top_lat, left_top_lon, right_bottom_lat, right_bottom_lon, num_features, output_shapefile, feature_type, atr_count, num_points):
    if num_points < 4:
        raise ValueError("num_points must be at least 4")

    # UTF-8エンコードを設定する
    os.environ['SHAPE_ENCODING'] = 'UTF-8'
    interval = 300

    # num_featuresを作成するための行と列の数を計算する
    grid_size = math.ceil(math.sqrt(num_features))
    rows = grid_size
    cols = grid_size

    # 各マス目またはステップの大きさを計算してポイントを獲得する
    lat_step = (left_top_lat - right_bottom_lat) / rows
    lon_step = (right_bottom_lon - left_top_lon) / cols

    # Shapefileを作成する
    driver = ogr.GetDriverByName("ESRI Shapefile")
    ds = driver.CreateDataSource(output_shapefile)
    geom_type = ogr.wkbPolygon if feature_type == 'polygon' else ogr.wkbPoint
    layer = ds.CreateLayer("features", geom_type=geom_type)

    # 座標参照系(CSR)を定義する WGS84
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromEPSG(4326)  # WGS84

    # 座標参照系(CSR)を .prj に保存する
    prj_path = os.path.splitext(output_shapefile)[0] + ".prj"
    with open(prj_path, 'w') as prj_file:
        prj_file.write(spatial_ref.ExportToWkt())

    for i in range(1, atr_count):
        if i % interval == 0:
            key = f"Atr{i:07}"
            layer.CreateField(ogr.FieldDefn(key, ogr.OFTInteger))
        else:
            key = f"Atr{i:07}"
            field_defn = ogr.FieldDefn(key, ogr.OFTString)
            field_defn.SetWidth(254)
            layer.CreateField(field_defn)

    # 最大文字れるの候補
    str_types = [generate_random_string(254), generate_random_string(254), generate_random_string(254), generate_random_string(254)]

    created_features = 0
    for row in range(rows):
        for col in range(cols):
            if created_features >= num_features:
                break

            if feature_type == 'polygon':
                ring = ogr.Geometry(ogr.wkbLinearRing)

                # 正方形の辺の点を計算する
                points = []
                lat1, lon1 = left_top_lat - row * lat_step, left_top_lon + col * lon_step
                lat2, lon2 = left_top_lat - row * lat_step, left_top_lon + (col + 1) * lon_step
                lat3, lon3 = left_top_lat - (row + 1) * lat_step, left_top_lon + (col + 1) * lon_step
                lat4, lon4 = left_top_lat - (row + 1) * lat_step, left_top_lon + col * lon_step

                # 上端
                for i in range(num_points // 4):
                    lat = lat1
                    lon = lon1 + i * (lon2 - lon1) / (num_points // 4)
                    points.append((lon, lat))

                # 右端
                for i in range(num_points // 4):
                    lat = lat2 + i * (lat3 - lat2) / (num_points // 4)
                    lon = lon2
                    points.append((lon, lat))

                # 下端
                for i in range(num_points // 4):
                    lat = lat3
                    lon = lon3 - i * (lon3 - lon4) / (num_points // 4)
                    points.append((lon, lat))

                # 左端
                for i in range(num_points // 4):
                    lat = lat4 - i * (lat4 - lat1) / (num_points // 4)
                    lon = lon4
                    points.append((lon, lat))

                # ポイントを追加
                for lon, lat in points:
                    ring.AddPoint(lon, lat)

                ring.CloseRings()

                polygon = ogr.Geometry(ogr.wkbPolygon)
                polygon.AddGeometry(ring)

                feature = ogr.Feature(layer.GetLayerDefn())
                feature.SetGeometry(polygon)
            else:
                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(left_top_lon + (col + 0.5) * lon_step, left_top_lat - (row + 0.5) * lat_step)

                feature = ogr.Feature(layer.GetLayerDefn())
                feature.SetGeometry(point)

            for i in range(1, atr_count):
                if i % interval == 0:
                    key = f"Atr{i:07}"
                    feature.SetField(key, random.randint(1, 100))
                else:
                    key = f"Atr{i:07}"
                    feature.SetField(key, random.choice(str_types))

            layer.CreateFeature(feature)
            feature = None

            created_features += 1

    ds = None

    cpg_path = os.path.splitext(output_shapefile)[0] + ".cpg"
    with open(cpg_path, 'w', encoding='utf-8') as cpg_file:
        cpg_file.write('UTF-8')

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Create Shapefile with specified features.")
    parser.add_argument("left_top_lat", type=float, help="Left top latitude")
    parser.add_argument("left_top_lon", type=float, help="Left top longitude")
    parser.add_argument("right_bottom_lat", type=float, help="Right bottom latitude")
    parser.add_argument("right_bottom_lon", type=float, help="Right bottom longitude")
    parser.add_argument("num_features", type=int, help="Number of features to create")
    parser.add_argument("output_shapefile", type=str, help="Output Shapefile path")
    parser.add_argument("feature_type", choices=['polygon', 'point'], help="Type of features to create (polygon or point)")
    parser.add_argument("--atr_count", type=int, help="Number of attributes(Max 256)", default=30)
    parser.add_argument("--num_points", type=int, help="Number of vertices in a Polygon", default=4)

    args = parser.parse_args()

    create_shapefile(args.left_top_lat, args.left_top_lon, args.right_bottom_lat, args.right_bottom_lon, args.num_features, args.output_shapefile, args.feature_type, args.atr_count, args.num_points)