On python, how to generate an array of coordinates between two points

Question

Perhaps I'm looking to a solved question, but while searching I didn't find any solution. I need to build a matrix of lat, lon values, based on two points, and a sampling rate (distance between points).

A possibility is implement a solution using numpy, with all the rules. However some exception treatment is needed, specially if my region is on an pole or mid-pacific zones, where a WGS-84 coordinate based matrix would change from 179.9 to -179.9, for example.

I would like to know if already exists a library that is able to generate interactively an array within those two points, and able to control those issues.

Answer 1

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

def pointiterator(fra,til,steps):
    '''
    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
        else:
            yield val
        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
xx=np.fromiter(xiter,dtype=np.float)
yy=np.fromiter(yiter,dtype=np.float)
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
print "x,y"
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.