Certainly doable in numpy. There may well be a better way of doing this, but I put the wrap-around logic in a generator method and used this to build a matrix of coordinates. This is then (un)ravelled into a series of x,y coordinates. Those can be saved to csv and loaded in as a delimited layer.
This code uses a step count, rather than a step size, but you should be able to change the pointiterator() method to suit.
import numpy as np
generator, like range() but uses number of steps,
and handles edge cases (like -180 to +180 crossover)
val = fra
if til < fra:
til += 360.0
stepsize = (til - fra)/steps
while val < til + stepsize:
if (val > 180.0):
yield val - 360.0
val += stepsize
# edge case - crosses antimeridian
xiter = pointiterator(150.0,-120.0,10)
yiter = pointiterator(-20.0,60.0,10)
# standard case
xiter = pointiterator(-20.0,70.0,10)
yiter = pointiterator(-10.0,20.0,10)
# make grid
xo, yo = np.meshgrid(xx,yy,indexing='xy')
# make array of tuples
tups = np.rec.fromarrays([xo,yo], names='x,y')
# dump points, load into qgis as delimited
for x,y in np.ravel(tups):
print "%f,%f" % (x,y)
Taking that into qgis as a delimited layer and showing on a web map...
That shows a standard case, and the edge case where the points overlap the antimeridian.
It doesn't handle the edge case at the poles. 91 degrees north doesn't wrap round to 89 degrees south ;) But it will work across the antimeridian. For the polar case you might need to tweak the generator to clip out any values > 90 or < -90.