0

I have this toy example that uses numpy and an external package called PyProj. Lat and lon are 2D arrays that contain coordinates of some domain specific information. What I want to do is to calculate distance on a sphere from a central point that I choose arbitrarily. The shape of lat_0 and lon_0 is

(2000,1)

But the API call inv does not like that. I get a run time error -

 RuntimeError: Buffer lengths not the same

It wants an array of shape 

(2000,50).

So I want lat_0 and lon_0 to be the same shape as lon and lat with all constant values which is the central latitude and longitude. What is the most efficient way to increase the columns of lon_0 and lat_0 and filling it with the central value so that it is the same shape as lon and lat without using for loops? 

import numpy as np
from pyproj import Geod

lat = np.empty((2000,50))
lat.fill(1)
lon = np.empty((2000,50))
lon.fill(1)


center = int(np.floor(len(lon[-1]) / 2.))
lon_0 = lon[:,center][...,np.newaxis]
lat_0 = lat[:,center][...,np.newaxis]


g = Geod(ellps='WGS84')

distance = g.inv(lon,lat,lon_0,lat_0,radians=True)
gansub
  • 1,164
  • 3
  • 20
  • 47

2 Answers2

2

The most efficient way would probably be np.broadcast_arrays. This creates views of the smaller arrays without enlarging the data buffer. Generic example:

    >>> A = np.arange(10).reshape(2, 5)
>>> A
array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])
>>> B = np.c_[:2]
>>> B
array([[0],
       [1]])
>>> C = np.arange(5)
>>> D = 7

>>> np.broadcast_arrays(A, B)
[array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]]), array([[0, 0, 0, 0, 0],
       [1, 1, 1, 1, 1]])]
>>> np.broadcast_arrays(A, C)
[array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]]), array([[0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4]])]
>>> np.broadcast_arrays(A, D)
[array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]]), array([[7, 7, 7, 7, 7],
       [7, 7, 7, 7, 7]])]

To see that the data is shared:

>>> AA, BB = np.broadcast_arrays(A, B)
>>> BB
array([[0, 0, 0, 0, 0],
       [1, 1, 1, 1, 1]])
>>> BB[0,0] = 3
>>> BB
array([[3, 3, 3, 3, 3],
       [1, 1, 1, 1, 1]])
>>> B
array([[3],
       [1]])
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • yours maybe the more efficient answer but I am a newbie to numpy broadcast_arrays. Can you provide some hints on the relation between your A, B and my lon,lat and lon_0 and lat_0 so that I can come up with a toy example ? – gansub Apr 11 '18 at 10:45
  • 1
    @gansub The shape of `B` is (2, 1), so that would correspond to your `lat_0`, `lon_0`; the shape of `A` is (2, 5), that would correspond to your `lat` and `lon`. If you want to deeply understand `broadcast_arrays` have a look at for example `BB.strides`. – Paul Panzer Apr 11 '18 at 11:45
1

Not 100% sure I understood what you need, my suggestion seems uneficient but to copy an array multiple times along a specified axis, you could use numpy repeat

So in your case you could do 

lon_0 = np.repeat(lon_0, 50, axis=1)
Christian
  • 1,162
  • 12
  • 21