5

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.

  • do you want to sample points along a line between two coordinates? Or do the coordinates mark the diagonal corners of a rectangular region, and you want a regular grid with points sampled in both axes? – Steven Kay Oct 22 '15 at 19:30
  • It would be the second option. The idea is generate a matrix of centroids. – Ricardo Barros Lourenço Oct 22 '15 at 20:09
5

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...

enter image description here

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.

  • Thanks Steven :) I think this addresses a lot of my issues. For pole, I would do some sort of domain validation, as you suggested. – Ricardo Barros Lourenço Oct 22 '15 at 20:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.