大量の地物数のShapefileを作成する機会があったので作成したプログラムをご紹介します。
サイズの小さいものや大きなもの地物数等いろいろな条件でShapefileが作成できるので、テスト用にご活用ください。
Shapefileは使用するソフトウェアや環境によって扱えるファイルサイズが異なる場合がありますが、一般的に使われている2GBを超えるShapefileも作成できることを確認しています。
プログラムの仕様
・Shapefileを作成する範囲を緯度経度で指定できる
・PolygonまたはPointのShapefileどちらをつくるか指定できる
・作成する地物数を指定できる
・作成する属性数を指定できる
・Polygonの場合1Polygonの頂点数を指定できる
・文字コードはUTF-8で作成する
プログラムの使い方
プログラムに指定する引数は下記の通りです。
引数1:作成する範囲左上の緯度
引数2:作成する範囲左上の経度
引数3:作成する範囲右下の緯度
引数4:作成する範囲右下の経度
引数5:作成する地物数
引数6:出力ファイルパス
引数7:作成する地物の種類(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)