如何从python中的中心点创建坐标网格？

Question

有没有一种方法可以从中心坐标中提取坐标网格（紫色点），例如每个坐标之间的距离为 100 米？

例如输入经纬度中心坐标（红点）：

lat = 40.8264859
lon = -3.6805897

在一个函数中，然后返回一个列表或数组的列表，其中包含相应的 lat、lon 网格，彼此相距 100 米，以及距中心点的最大坐标数（在此图像情况下，该值将是2)。例如：

def get_grid_of_coordinates_from_center(lat, lon, meters_between_coor, coors_away_from_center):
    ...
    return list_of_lists

Answer 1

试试这种矢量化方法。用米抵消 lat-long 的方程式灵感来自 this link on stack exchange -

import numpy as np

lat, lon = 50, -10 #center coordinate
dist, coors = 100, 2 #meters, num coordinates in each direction

#Creating the offset grid
mini, maxi = -dist*coors, dist*coors
n_coord = coors*2+1
axis = np.linspace(mini, maxi, n_coord)
X, Y = np.meshgrid(axis, axis)


#avation formulate for offsetting the latlong by offset matrices
R = 6378137 #earth's radius
dLat = X/R
dLon = Y/(R*np.cos(np.pi*lat/180))
latO = lat + dLat * 180/np.pi
lonO = lon + dLon * 180/np.pi

#stack x and y latlongs and get (lat,long) format
output = np.stack([latO, lonO]).transpose(1,2,0)
output.shape

(5,5,2)

让我们绘图以查看点的网格并确认这些点正确分布在中心 lat long 周围。

import matplotlib.pyplot as plt

points = output.reshape(-1,2)
x = points[:,0]
y = points[:,1]

plt.scatter(x,y)              #<- plot all points
plt.scatter(50,-10,color='r') #<- plot the center lat long

---

解释

为了便于理解，我将上述步骤分解为单独的函数（并使用np.vectorize）。

1. 首先，您需要一个矢量化函数，它接收（纬度、经度）并在 X 和 Y 方向上偏移米。
2. 您需要创建一个offset_grid
   • 您可以首先使用np.linspace 创建一个范围，从头到尾具有所需的点数
   • 接下来，您可以使用它来创建一个带有 X、Y 矩阵的网格，使用 np.meshgrid
3. 使用矢量化函数和offset_grid，您可以对每个值应用经纬度偏移，然后重新整形以获得 (n,n) 预期矩阵中每个点的 (x,y) 值。

让我们从一个用一些 x 和 y 米偏移 lat、long 的函数开始。

#Offset any lat long by x, y meters
def lat_long_offset(lat, lon, x, y):
    '''
    lat, lon : Provide lat lon coordinates
    x, y : Provide offset of x and y on lat and long respectively
           This needs to be in meters!
           
    The approximation is taken from an aviation formula from this stack exchange 
    https://gis.stackexchange.com/questions/2951/algorithm-for-offsetting-a-latitude-longitude-by-some-amount-of-meters
    '''
       
    #Earth’s radius, sphere
    R=6378137

    #Coordinate offsets in radians
    dLat = x/R
    dLon = y/(R*np.cos(np.pi*lat/180))

    #OffsetPosition, decimal degrees
    latO = lat + dLat * 180/np.pi
    lonO = lon + dLon * 180/np.pi
    
    return latO, lonO

#Create a vectorized offset function
lat_long_offset_vec = np.vectorize(lat_long_offset)

一旦我们准备好这个函数，我们就可以简单地在偏移网格上工作，我们需要在所有方向上应用偏移，以获得经纬度格式的相关点坐标。

#Create offset_grid and return coordinates
def get_mesh(lat, lon, dist, coors):
    #calculate min and max range for coordinates over an axis
    mini, maxi = -dist*coors, dist*coors
    
    #calculate number of points over an axis
    n_coord = coors*2+1
    
    #create an axis from min to max value with required number of coordinates
    axis = np.linspace(mini, maxi, n_coord)
    
    #create an "offset_grid" for X and Y values for both axis. 
    X, Y = np.meshgrid(axis, axis)

    #calcualte offset coordinates for "offset_grid" in meters
    mesh = lat_long_offset_vec(lat, lon, X, Y)
    
    #Transpose to get the (x,y) values for the offset_grid's shape
    mesh_x_y_format = np.stack(mesh).transpose(1,2,0)
    return mesh_x_y_format

output = get_mesh(50, -10, 100, 2)

print('Shape of output grid:', output.shape)
print('Note: 2 values (x,y) for each point in the expected (5,5) grid')

print('')
print('Output coordinates')
print(output)

Shape of output grid: (5, 5, 2)
Note: 2 values (x,y) for each point in the expected (5,5) grid

Output coordinates
[[[ 49.99820337 -10.00279506]
  [ 49.99910168 -10.00279506]
  [ 50.         -10.00279506]
  [ 50.00089832 -10.00279506]
  [ 50.00179663 -10.00279506]]

 [[ 49.99820337 -10.00139753]
  [ 49.99910168 -10.00139753]
  [ 50.         -10.00139753]
  [ 50.00089832 -10.00139753]
  [ 50.00179663 -10.00139753]]

 [[ 49.99820337 -10.        ]
  [ 49.99910168 -10.        ]
  [ 50.         -10.        ]  #<- Notice the original coordinate at center
  [ 50.00089832 -10.        ]
  [ 50.00179663 -10.        ]]

 [[ 49.99820337  -9.99860247]
  [ 49.99910168  -9.99860247]
  [ 50.          -9.99860247]
  [ 50.00089832  -9.99860247]
  [ 50.00179663  -9.99860247]]

 [[ 49.99820337  -9.99720494]
  [ 49.99910168  -9.99720494]
  [ 50.          -9.99720494]
  [ 50.00089832  -9.99720494]
  [ 50.00179663  -9.99720494]]]

Answer 2

可能需要检查一些极端情况（我使用了易于使用的数字），但这似乎可行

import matplotlib.pyplot as plt
import numpy as np

mid_x, mid_y, max_x, max_y, min_x, min_y,step_size,num_points = 5, 5, 11, 11, 0, 0, 1, 10
x_range = np.concatenate((np.arange(min_x, mid_x, step_size), np.arange(mid_x, max_x, step_size)))
y_range = np.concatenate((np.arange(min_y, mid_y, step_size), np.arange(mid_y, max_y, step_size)))
data = np.array([[x, y] for x in x_range for y in y_range])  # cartesian prod
plt.scatter(data[:, 0], data[:, 1])
plt.scatter(mid_x, mid_y, color='r')
plt.show()

情节：