0

I have a Numpy array of shape (30000, 2000, 2) where each row is supposed to have 2000 latitude and longitude points with equal distance from each other and in the same direction. The problem is, each row has missing points. For each row, just the latitude and longitude of the start points and the endpoints (in CRS WSG84) are known and the rest of the points are missing. For example a single row can be like this:

latlonPoints[0,:] = [(48.778767, -123.903275),?,?,?,...,(49.887672, -122.49999)]

Using the start and endpoints for each row, we can easily calculate the distance (sometimes a few kilometers) and also the angle (bearing) between the two points. Then using the distance and bearing we can generate the rest of the points using geopy or other methods. But, iteration using for loop is too slow for all the rows and all the elements for every single row. What will be the fastest way in python to calculate the missing points? Here is what I have implemented and it works fine but slow:

    for idx, row in enumerate(latlonPoints):
    lat_start = row[0]
    long_start = row[1]
    coords1 = (lat_start, long_start)
    lat_end = row[2]
    long_end = row[3]
    bearing, back_azimuth, distance = geodesic1.inv(long_start, lat_start, long_end, lat_end)
    step_size = distance/(latlonPoints.shape[1] -1)
    latlon_row = []
    for i in range (0,latlonPoints.shape[1]):
        dist = i * step_size
        destination = geodesic(kilometers=dist / 1000).destination(coords1, bearing)
        latlon_row.append((destination.latitude,destination.longitude))
    latlonPoints[idx,:] = np.array(latlon_row)
Cheeta_2019
  • 1
  • 3
  • 1
    list comprehension, `numpy` loop or something or a mathematical formula (the fastest), otherwise use some faster language – Matiiss Oct 18 '21 at 23:02
  • 1
    [This](https://gis.stackexchange.com/questions/311362/interpolating-orthodrome-between-two-lon-lat-points-in-python) answer can be helpful. The library seems to be able to interpolate 2000 points in under one millisecond on my machine. – hilberts_drinking_problem Oct 18 '21 at 23:26
  • If you're really just spanning a few kilometers, there's no need to use spherical coordinates. Just do linear interpolation. The difference will not be noticeable. – Tim Roberts Oct 19 '21 at 03:08
  • @hilberts_drinking_problem "pyproj.Geod(ellps='WGS84')" does the trick quite fast. Thanks – Cheeta_2019 Oct 19 '21 at 22:19

1 Answers1

0

As @hilberts_drinking_problem answered the question using poproj and Geod, I can generate multiple points Given a single initial point and terminus point. It returns a list of longitude/latitude pairs describing npts equally spaced intermediate points along the geodesic between the initial and terminus points.

from pyproj import Geod

extra_points = pyproj.Geod(ellps='WGS84')
Cheeta_2019
  • 1
  • 3