If you are serious about sampling uniformly within a distance of a
particular point then you need account for the curvature of the
ellipsoid. This is pretty easy to do using a rejection technique. The
following Matlab code implements this.
function [lat, lon] = geosample(lat0, lon0, r0, n)
% [lat, lon] = geosample(lat0, lon0, r0, n)
%
% Return n points on the WGS84 ellipsoid within a distance r0 of
% (lat0,lon0) and uniformly distributed on the surface. The returned
% lat and lon are n x 1 vectors.
%
% Requires Matlab package
% http://www.mathworks.com/matlabcentral/fileexchange/39108
todo = true(n,1); lat = zeros(n,1); lon = lat;
while any(todo)
n1 = sum(todo);
r = r0 * max(rand(n1,2), [], 2); % r = r0*sqrt(U) using cheap sqrt
azi = 180 * (2 * rand(n1,1) - 1); % sample azi uniformly
[lat(todo), lon(todo), ~, ~, m, ~, ~, sig] = ...
geodreckon(lat0, lon0, r, azi);
% Only count points with sig <= 180 (otherwise it's not a shortest
% path). Also because of the curvature of the ellipsoid, large r
% are sampled too frequently, by a factor r/m. This following
% accounts for this...
todo(todo) = ~(sig <= 180 & r .* rand(n1,1) <= m);
end
end
This code samples uniformly within a circle on the azimuthal equidistant
projection centered at lat0, lon0. The radial, resp. azimuthal,
scale for this projection is 1, resp. r/m. Hence the areal
distortion is r/m and this is accounted for by accepting such points
with a probability m/r.
This code also accounts for the situation where r0 is about half the
circumference of the earth and avoids double sampling nearly antipodal
points.
You should try the 'RandomPointsInPolygon(geom geometry, num_points integer)` function.
Add an additional geometry column, and update it using RandomPointsInPolygon. The first parameter is your existing geometry, and the second is your PNTSneeded column.
There are two question on SO related with RandomPointsInPolygon:
How to create random points in a polygon in PostGIS
Postgis random point inside a polygon
Start with a smaller sample, if possible, using a WHERE clause in the UPDATE statement.
The actual code would be:
-- copy and paste the RandomPointsInPolygon function
CREATE OR REPLACE FUNCTION RandomPointsInPolygon(geom geometry, num_points integer)
RETURNS SETOF geometry AS
$BODY$DECLARE
(...)
-- if you use EPSG:4326
SELECT AddGeometryColumn('polygons', 'points', 4326, 'MULTIPOINT', 2);
-- update your table polygons, using existing geom and PNTSneeded
UPDATE polygons
SET points = (SELECT ST_Union(manypoints)
FROM RandomPointsInPolygon("geom", "PNTSneeded") AS manypoints);
Best Answer
I made and tested the functions I need. For point generation I used a simplified version of the functions in trac.osgeo.org/postgis/wiki/UserWikiRandomPoint:
I don't need the position to be extra accurate, so I use approximate coordinates for Czech Republic.
I don't need to generate lines now. For the polygon I started from scratch, using the algorithm I described in the question (building the polygon around a quasicentroid):
I picked few of the resulting geometries for the following image:
I don't mind when geometries intersect each other; otherwise there would be a check for this too.