Generate zipped shapefiles of random points

random_points.py

#!/usr/bin/env python3
"""Generate zipped shapefiles of random points inside a square area.

Only dependency: GDAL python bindings (osgeo).

Given a center (lon/lat EPSG:4326), the side of a square in km, the number
of points per shapefile and the number of shapefiles, it creates:

    ./points_<npoints>_<lon>_<lat>_<crs>/
        <npoints>_<crs>_1.zip
        <npoints>_<crs>_2.zip
        ...
Each zip contains the companion files of one shapefile (.shp .shx .dbf .prj).

Copyright 2026 Stefano Cudini https://github.com/stefanocudini
"""

import argparse
import os
import random
import string
import sys
import zipfile

from osgeo import ogr, osr

ogr.UseExceptions()
osr.UseExceptions()

EPSG_WGS84 = 4326

SHP_EXTS = (".shp", ".shx", ".dbf", ".prj", ".cpg")


def utm_epsg_for(lon, lat):
    """EPSG code of the UTM zone containing lon/lat."""
    zone = int((lon + 180) // 6) + 1
    return (32600 if lat >= 0 else 32700) + zone


def build_transforms(lon, lat):
    """Return (to_utm, to_wgs84) coordinate transformations."""
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(EPSG_WGS84)
    wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    utm = osr.SpatialReference()
    utm.ImportFromEPSG(utm_epsg_for(lon, lat))
    utm.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    return (osr.CoordinateTransformation(wgs84, utm),
            osr.CoordinateTransformation(utm, wgs84))


def random_points(lon, lat, side_km, n):
    """n random (lon, lat) points uniform in a square of side_km centered on lon/lat."""
    to_utm, to_wgs84 = build_transforms(lon, lat)
    cx, cy, _ = to_utm.TransformPoint(lon, lat)
    half = side_km * 1000.0 / 2.0
    pts = []
    for _ in range(n):
        x = random.uniform(cx - half, cx + half)
        y = random.uniform(cy - half, cy + half)
        plon, plat, _ = to_wgs84.TransformPoint(x, y)
        pts.append((plon, plat))
    return pts


def write_shapefile(shp_path, points):
    """Write points (lon/lat) to an ESRI Shapefile in EPSG:4326."""
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(EPSG_WGS84)
    srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    driver = ogr.GetDriverByName("ESRI Shapefile")
    ds = driver.CreateDataSource(shp_path)
    layer = ds.CreateLayer(os.path.splitext(os.path.basename(shp_path))[0],
                           srs, ogr.wkbPoint)
    layer.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
    name_field = ogr.FieldDefn("name", ogr.OFTString)
    name_field.SetWidth(10)
    layer.CreateField(name_field)

    for i, (plon, plat) in enumerate(points):
        feature = ogr.Feature(layer.GetLayerDefn())
        feature.SetField("id", i)
        feature.SetField("name", "".join(random.choices(string.ascii_lowercase, k=10)))
        geom = ogr.Geometry(ogr.wkbPoint)
        geom.AddPoint(plon, plat)
        feature.SetGeometry(geom)
        layer.CreateFeature(feature)
        feature = None

    ds = None  # flush to disk


def zip_shapefile(shp_path, zip_path):
    """Zip all companion files of shp_path into zip_path and remove them."""
    base = os.path.splitext(shp_path)[0]
    with zipfile.ZipFile(zip_path, "w", zipfile.ZIP_DEFLATED) as zf:
        for ext in SHP_EXTS:
            part = base + ext
            if os.path.exists(part):
                zf.write(part, os.path.basename(part))
    for ext in SHP_EXTS:
        part = base + ext
        if os.path.exists(part):
            os.remove(part)


def lonlat(value):
    """Parse a "lon,lat" string into a (lon, lat) float tuple."""
    try:
        lon, lat = (float(v) for v in value.split(","))
    except ValueError:
        raise argparse.ArgumentTypeError(
            '"{}" is not a valid "lon,lat" pair'.format(value))
    return lon, lat


def main():
    parser = argparse.ArgumentParser(
        description="Generate zipped shapefiles of random points in a square area.")
    parser.add_argument("--center", type=lonlat, default=(0.0, 0.0),
                        help='center as "lon,lat" (EPSG:4326), default "0,0"')
    parser.add_argument("--side-km", type=float, default=1000.0,
                        help="side of the square in km, default 1000")
    parser.add_argument("--points", type=int, default=1000,
                        help="number of points per shapefile, default 1000")
    parser.add_argument("--files", type=int, default=1,
                        help="number of shapefiles to generate, default 1")
    args = parser.parse_args()

    lon, lat = args.center
    out_dir = "./points_{n}_{lon}_{lat}_{crs}".format(
        n=args.points, lon=lon, lat=lat, crs=EPSG_WGS84)
    os.makedirs(out_dir, exist_ok=True)

    for i in range(1, args.files + 1):
        name = "random_points{n}_{crs}_{i}".format(n=args.points, crs=EPSG_WGS84, i=i)
        shp_path = os.path.join(out_dir, name + ".shp")
        zip_path = os.path.join(out_dir, name + ".zip")

        pts = random_points(lon, lat, args.side_km, args.points)
        write_shapefile(shp_path, pts)
        zip_shapefile(shp_path, zip_path)
        print(zip_path)

    print("Done: {} shapefile zip(s) in {}".format(args.files, out_dir),
          file=sys.stderr)


if __name__ == "__main__":
    main()