python - 如何从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

试试这种矢量化方法。用米抵消纬度经度的方程的灵感来自 this link在堆栈交换上 -

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)

让我们绘制以查看点网格并确认点是否正确分布在中心纬度周围。

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 )。

首先，您需要一个矢量化函数，它接收 (lat, long) 并在 X 和 Y 方向上以米为单位偏移它。

您需要创建一个 offset_grid

您可以从 np.linspace 开始使用所需的点数创建从开始到结束的范围

接下来，您可以使用 np.meshgrid 创建一个包含 X、Y 矩阵的网格。

使用矢量化函数和 offset_grid您可以对每个值应用经纬度偏移，然后 reshape 以获得 (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]]]