Hướng dẫn python azimuth - phương vị trăn

Tôi có hai tọa độ trên thế giới và tôi cần tính toán góc phương vị giữa chúng, có một số chức năng trong Python3 để tính toán nó cho tôi chỉ bằng cách nhập hai tọa độ không?Azimuth between them, Is there some function in Python3 to calculate it for me just by entering the two coordinates?

Các tọa độ là X, Y, Z tính bằng mét và trong hệ tọa độ "UTM WGS84". Tôi biết có một số cách dài để làm điều đó, tôi chỉ không biết làm thế nào và tôi thích một cách ngắn để làm điều đó. Cảm ơn.UTM WGS84" Coordinate System. I know there is some long way to do it I just dont know how and I prefer a short way to do it. Thanks.

Ví dụ về hai tọa độ (tính bằng mét):

X1 = 10.3, y1 = 20.4, z1 = 4,8

X2 = 14.2, y2 = 35.1, z2 = 5.0

Nếu có thêm một số thông tin tôi cần thêm PLZ cho tôi biết và tôi sẽ thêm nó.

Sau vài ngày tìm kiếm và thử nghiệm, cuối cùng tôi cũng có được điều này để làm việc. Câu hỏi đầu tiên là liệu có thể tính toán Azimuth bằng tọa độ x, y pixel của hình ảnh hay không. Câu trả lời là có. Theo câu trả lời ở đây, bạn có thể tính toán góc phương vị/ổ trục trong hình ảnh 2D. Đây là cách mã trông trong Python:
The first question was whether it was possible to calculate azimuth using an image's x,y pixel coordinate.
The answer is Yes.
Following the answer here, you can calculate the azimuth/bearing in a 2D image.
This is how the code looks in python:

def calculate_bearing(ax, ay, bx, by):
    """Computes the bearing in degrees from the point A(a1,a2) to
    the point B(b1,b2). Note that A and B are given in terms of
    screen coordinates.

    Args:
        ax (int): The x-coordinate of the first point defining a line.
        ay (int): The y-coordinate of the first point defining a line.
        bx (int): The x-coordinate of the second point defining a line.
        by (int): The y-coordinate of the second point defining a line.

    Returns:
        float: bearing in degrees
    """
    
    TWO_PI = math.pi * 2
    # if (a1 = b1 and a2 = b2) throw an error 
    theta = math.atan2(bx - ax, ay - by)
    if (theta < 0.0):
        theta += TWO_PI
    return math.degrees(theta) 

Ngoài ra, có ý chính này, nơi nó tính toán hướng của một dòng, theo cách tương tự, được đề xuất ở đây.

# https://gist.github.com/aarondandy/604912
# for calculating the bearing/azimuth from image pixel coordinates
def calculate_line_direction (ax, ay, bx, by):
    """Calculates the north relative clockwise orientation of a line.

    Args:
        ax (int): The x-coordinate of the first point defining a line.
        ay (int): The y-coordinate of the first point defining a line.
        bx (int): The x-coordinate of the second point defining a line.
        by (int): The y-coordinate of the second point defining a line.

    Returns:
        float: The direction of the line in degree or nan on failure.
    """
    dx = bx - ax
    dy = by - ay
    m = math.sqrt((dx * dx) + (dy * dy))
    degree = math.degrees(math.atan2(dx / m, dy / m))
    return degree if m>0 else math.nan

Ngoài ra, tôi phát hiện ra rằng việc tính toán ổ trục trong trường hợp của tôi có thể được thực hiện theo một cách tương đối đơn giản.

Vì góc phương vị là tiêu đề được đo theo chiều kim đồng hồ từ phía bắc và hình ảnh của chúng tôi đã được căn chỉnh về phía bắc, điểm bên phải sẽ có góc 90 độ (theo chiều kim đồng hồ) và điểm dưới đây sẽ có góc 180 độ. Chúng ta có thể mã hóa những thứ này hoặc tính toán chúng giống như phương pháp tôi vừa đề cập. Đây là một con số đơn giản mà tôi đã đưa ra để thể hiện điều này.
We can hardcode these or calculate them like the method I just mentioned. this is a simple figure I came up with to show this.

    #  
    # |  -    
    # |      -
    # |         -
    # |           -
    # |             -
    # |               -
    # |                 -
    # |                  -
    # |                   -
    # |                    -
    # |                     -
    # point-a(0)--------point-b (90)
    # |                     - 
    # |                    -
    # |                   -
    # |                  -
    # |                -
    # |              -
    # |            -
    # |         -
    # |      -
    # |   -
    # point-c(180)

Trước khi tôi tiếp tục, tôi muốn nói về tính toán tỷ lệ pixel. Không có tính toán tỷ lệ pixel chính xác, mọi thứ đều có trong tĩnh mạch. Tìm kiếm cách tính thang đo pixel kết thúc với rất nhiều kết quả không hoạt động! Giống như người được giải thích ở đây:

# This doesn't work! don't use it!
const TILE_SIZE = 256;          // Tile size of base tile.
const CIRCUMFERENCE = 40075016; // Circumference at the equator in meters.
const scale = 1 << zoomLevel;   // Bitshifting is easy shorthand for Math.pow(2, zoomLevel)
const metersPerPixel = CIRCUMFERENCE / (TILE_SIZE * scale);

Vấn đề là, nếu bạn sử dụng đoạn trích ở trên, bạn sẽ nhận được những sai lệch lớn dựa trên nơi bạn sử dụng trên thế giới để tính toán thang đo pixel. Về cơ bản khi nói đến phép chiếu Mercator, nó tạo ra một thang đo ngang khác nhau ở các vĩ độ khác nhau và do đó thang đo pixel khác nhau dựa trên vĩ độ. Điều dường như nhận được một xấp xỉ rất tốt là
Basically when it comes to Mercator projection, it produces a different horizontal scale at different latitudes and thus pixel scale differs based on the latitudes.
What seems to get a very decent approximation is

# https://groups.google.com/g/google-maps-js-api-v3/c/hDRO4oHVSeM
# https://gis.stackexchange.com/a/127949/191141
def calculate_pixel_scale(lat, zoom):
    """Calculates and returns the pixel scale in meters at the specified zoom level.
        This function shows what a pixel represents in real world. using this you can get
        the real world estimate of distance traversed, or object sizes, etc in meters. 
        This is based on the GoogleMap's specifications (zoom, tile size, etc). 

    Args:
        lat (float): the latitude at which you want to calculate the pixel scale
        zoom (int): the zoom level at which you want to calculate the pixel scale

    Returns:
        float: pixel scale
    """
    scale = 156543.03392 * math.cos(lat * math.pi / 180) / math.pow(2, zoom)
    return scale

Bây giờ, bằng cách nào đó để có được điều này và tính toán vĩ độ/kinh độ tiếp theo, tôi cần phải có vĩ độ của vị trí đó! Có vẻ như một vòng tròn bây giờ! Để phá vỡ chu kỳ này, chúng ta có thể sử dụng tọa độ Google Ngói. Chúng tôi biết rằng tại các zoom khác nhau, chúng tôi có số lượng gạch khác nhau tạo thành bản đồ của chúng tôi. Những viên gạch này đều có kích thước 256x256, vì vậy nếu chúng ta biết một gạch duy nhất ở mức zoom được chỉ định trong bản đồ của chúng tôi, chúng ta có thể tính toán bất kỳ tọa độ gạch nào khác trong bản đồ của chúng tôi (vì chúng là 256x256, đó là một toán học đơn giản) vì vậy nếu chúng ta có thể nhận được Tọa độ gạch ở tọa độ pixel của bản đồ (0,0) của chúng tôi, chúng tôi vẫn ổn. Có nhiều cách để làm điều này, nhưng một cách đơn giản sẽ là tính toán tọa độ gạch từ một vĩ độ, kinh độ nhất định. Điều này có thể được thực hiện như thế này:
In order to break this cycle, we can use google tile coordinates. we know that at different zooms, we have different number of tiles that constitute our map.
These tiles are all 256x256 in size, so if we know a single tile at our specified zoom level in our map, we can calculate any other tile coordinates in our map (as they are 256x256, its a simple math)
So if we can get the tile coordinate at our map(0,0) pixel coordinate, we be fine. there are many ways to do this, but one simple way would be to calculate tile coordinate from a given latitude, longitude. This can be done like this:

# calculate tile coordinate based on lat/lon
# https://github.com/Coderx7/google-map-downloader/blob/master/downloader_1.1.py#L110
# note read the note!
def calculate_tile_coordinate(latitude, longitude, zoom):
    """Get google-style tile cooridinate from geographical coordinate
    Note: 
    Note that this will return a tile coordinate. and a tile coordinate
    is located at (0,0) or origin of a tile. 
    all tiles are 256x256. there are as many gps locations inside of a tile 
    that when used with this function will return the same exact tile coordinate
    (as they all reside in that tile obviously) in order to get the displacement 
    of your gps location, after calculating the tile coordinate, calculate tile's 
    (0,0) gps location from the newly achieved tile coordinate. then you have gps location
    at (0,0) and now can calculate howfar your initial gps location is from the 
    origin of the tile. 

    Args:
        latitude (float): Latitude
        longitude (float): Longittude
        zoom (int): zoom

    Raises:
        TypeError: [description]
        TypeError: [description]

    Returns:
        tuple(int,int): returns tile coordinate in the form of (y,x)
    """
    isnum = lambda x: isinstance(x, int) or isinstance(x, float)
    if not (isnum(longitude) and isnum(latitude)):
        raise TypeError("latitude and longitude must be int or float!")

    if not isinstance(zoom, int) or zoom < 0 or zoom > 22:
        raise TypeError("zoom must be int and between 0 to 22.")

    if longitude < 0:
        longitude = 180 + longitude
    else:
        longitude += 180
    longitude /= 360  # make longitude to (0,1)

    latitude = 85.0511287798 if latitude > 85.0511287798 else latitude
    latitude = -85.0511287798 if latitude < -85.0511287798 else latitude
    latitude = math.log(math.tan((90 + latitude) * math.pi / 360)) / (math.pi / 180)
    latitude /= 180  # make latitude to (-1,1)
    latitude = 1 - (latitude + 1) / 2  # make latitude to (0,1) and left top is 0-point

    num = 2 ** zoom
    y = math.floor(latitude * num)
    x = math.floor(longitude * num)
    return y,x

Bây giờ chúng ta có tọa độ gạch, chúng ta có thể kiểm tra các lớp pixel mong muốn của chúng ta cư trú, sau đó chúng ta có được tọa độ GPS của GP Vị trí GPS chính xác của tọa độ pixel được chỉ định. Ngoài ra, bây giờ chúng ta có thể có bất kỳ số gạch nào, bằng cách sử dụng gạch bên phải và gạch bên dưới gạch đích, chúng ta có thể tính toán tỷ lệ pixel chính xác cho x và trục y tương ứng (hiển thị bên dưới).
also now that we can have any tile number, by using the tile to the right and the tile below the target tile, we can calculate exact pixel scale for x, and y axis respectively(shown below).

Đây là mã nhanh và bẩn của tôi mà tôi đã chơi dí dỏm, nhưng trước đó, để trực quan hóa tốt hơn, điều này có thể giúp ích:

def get_distance(gps_coordinate1, gps_coordinate2):
    """Calculates the distance between to gps coordinates in meters.

    Args:
        coord1 (tuple(float, float)): latitude and longitude of point A
        coord2 (tuple(float, float)): latitude and longitude of point B

    Returns:
        float: distance is meters
    """
    return Geodesic.WGS84.Inverse(gps_coordinate1[0],
                                  gps_coordinate1[1], 
                                  gps_coordinate2[0], 
                                  gps_coordinate2[1])['s12']

def calculate_latitue_longitude_from_tile(row, column, zoom_level=18):
    """ 
        Given a tile coordinate, calculates the latitude and longitude at 
        the specified zoom level.

    Args:
        row (int): tile's row number. 
        column (int): tile's column number. 
        zoom_level (int, optional): the zoom at which level, map is collected. Defaults to 18.

    Returns:
        [type]: [description]
    """
    lon = column / math.pow(2.0, zoom_level) * 360.0 - 180
    n = math.pi - (2.0 * math.pi * row) / math.pow(2.0, zoom_level)
    lat = -((math.atan(math.sinh(n))) * -1) / math.pi * 180
    return lat, lon

def calculate_gps_from_pixel_offsets(x1=256, y1=0, offset_x=256, offset_y=0, target_tile=None, use_offsets_for_bearing=False, zoom=18):
    """
        calculates a new coordinate given an initial coordinate and pixel values x,y 
        In order to calculate a new coordinate from an existing one with a given distance 
        in meters, we need to calculate the bearing/azimuth. that is the direction at which
        we are headed needs to be known. for this we can calculate the angle between x,y 
        in our image. 

        the image we just talked about is our map. The map that was 
        from google earth. Knowing the tile numbers, we can calculate that locations coordinate
        since tiles are 256x256, and we created our map from these tiles, the (0,0) of our map
        corresponds to a tile number, we convert this tile coordinates to gps coorinate
        and we have our first coordinate ready. now we need to specify at what pixel location 
        we want to calculate the respective gps location.

        Suppose at (512,768), we calculate the angle between these pixels (row,col) with respec to
        north (googlemap is aligned to north) since we have coordinate at 0,0 of the map, other 
        locations will be on right(east) and down(south) direction. so calculating the azimuth from 
        something like (0,256) and (256,0) would apply to all other values of x,y in our map. 
        
        So we simply set x1,andy1 to 256 and 0 respectively. (0,256) will give us the azimuth for 
        x axis, and (256,0) will give us the azimuth for y axis, and using these two azimuths we 
        calculate two new locations/coordinates we use latitude from the one that uses y-azimuth 
        and use the longitude that was calculated using x-azimuth. this will give us the final 
        coordinate that we want. 

        as you can see, we also need the distance, offset_x and offset_y, they represent the 
        distance from our target tile pixel location. since they are in pixels
        we need to calculate the pixel scale at the zoom level we took the image, this will allow us
        to multiply the pixel offsets with this scale and get the metric distance. we can then
        use this metric distance to reach from the initial coordinate to the final one that we are after

    Args:
        x1 (int, optional): target x-pixel coord, we want to get its gps coord. Defaults to 256.
        y1 (int, optional): target y-pixel coord, we want to get its gps coord. Defaults to 0.
        offset_x (int, optional): pixel distance in x axis from tile's origin. Defaults to 256. it can be a value from 0-256
        offset_y (int, optional): pixel distance in y axis from tile's origin. Defaults to 0. it can be a value from 0-256
        target_tile (tuple(float,float), optional): the target tile coordinate that's used for calculating other coordinates for each pixel. Defaults to None.
        use_offsets_for_bearing (bool, optional): whether to use pixel coordinates for calculating azimuth instead of x1,y1. Defaults to False.
    """
    
    if use_offsets_for_bearing:
        x1, y1 = offset_x, offset_y

    # calculate azimuth/bearing for x,y individually
    azimuth_x = calculate_bearing(0, 0, x1, 0)
    azimuth_y = calculate_bearing(0, 0, 0, y1)

    # calculate the gps coordinate
    target_tile_gps = calculate_latitue_longitude_from_tile(target_tile[0], target_tile[1], zoom)

    # pixel scale at zoom 18
    # This calculates the azimuth, and is fairly accurate, however
    # the better way for calculating the azimuth is by using a separate
    # scaler for x,y respectively, this can be achieved by calculating the 
    # azimuth from 3 gps locations, one to the right and one down. 
    # these gps locations can be easily calculate by having one tile which
    # is our target tile, the tile that contains our pixel coordinate x,y
    # we then can simply calculate the other two tiles using simple math 
    # simply (tile_row, tile_col+1) and (tile_row+1, tile_col) will give us
    # the tiles, which we can then calculate the gps using google's formula
    # given before. this is done in our second method
    pixel_scale = calculate_pixel_scale(target_tile_gps[0], 18)

    # calculating the offset in meters based on pixel scale
    dsx_meters = offset_x * pixel_scale
    dsy_meters = offset_y * pixel_scale
    
    # if we set x1 and y1 to actual pixel locations(dsx,dsy), we need to remove the 180
    if not use_offsets_for_bearing:
        azimuth_x = 180 - azimuth_x 
        azimuth_y = 180 - azimuth_y 
   
    # to right/east direction
    g = Geodesic.WGS84.Direct(target_tile_gps[0], target_tile_gps[1], azimuth_x, dsx_meters)
    # to down/south direction
    g2 = Geodesic.WGS84.Direct(target_tile_gps[0], target_tile_gps[1], azimuth_y, dsy_meters)
    
    final_gps = g2['lat2'], g['lon2']

    print(f'\n---------------Method I-------------')
    print(f'distance offset in x axis: {dsx_meters:.2f}m.\ndistance offset in y axis: {dsy_meters:.2f}m.')
    print(f"GPS based on dsx: ({g['lat2']:.8f}, {g['lon2']:.8f}) az1: {azimuth_x:.2f}\N{DEGREE SIGN}")
    print(f"GPS based on dsy: ({g2['lat2']:.8f}, {g2['lon2']:.8f}) az2: {azimuth_y:.2f}\N{DEGREE SIGN}")
    print(f"GPS based on dsx-dsy: ({g2['lat2']:.8f}, {g['lon2']:.8f}).")
    print(f'Tile GPS: {target_tile_gps}')
    print(f"Final GPS:{final_gps}")
    print(f'Distance from tile gps to calculared gps location: {get_distance(target_tile_gps, final_gps):.2f}m.')

    return final_gps

def calculate_gps_from_pixel_offsets_accurate(target_tile, offset_x, offset_y, zoom=18):
   
    # calculate new tile's coordinate
    row, column = target_tile
    tile_right = (row, column+1)
    tile_down = (row+1, column)

    # calculate gps coordinates
    tile_gps = calculate_latitue_longitude_from_tile(row, column, zoom)
    tile_right_gps = calculate_latitue_longitude_from_tile(row, column+1, zoom)
    tile_down_gps = calculate_latitue_longitude_from_tile(row+1, column, zoom)

    # Define the path from 1 to 2
    lx = Geodesic.WGS84.InverseLine(tile_gps[0], tile_gps[1], tile_right_gps[0], tile_right_gps[1])
    ly = Geodesic.WGS84.InverseLine(tile_gps[0], tile_gps[1], tile_down_gps[0], tile_down_gps[1])
    
    # calculating pixel scale:
    # between tiles, there is a 256x256 distance in (x,y)
    pixel_scale_x = lx.s13/256
    pixel_scale_y = ly.s13/256

    # distance in meters in x axis
    dx_meters = offset_x*pixel_scale_x
    # distance in meters in y axis
    dy_meters = offset_y*pixel_scale_y
    
    # https://gis.stackexchange.com/questions/349458/
    direction_x = Geodesic.WGS84.Direct(tile_gps[0], tile_gps[1], lx.azi1, dx_meters)
    direction_y = Geodesic.WGS84.Direct(tile_gps[0], tile_gps[1], ly.azi1, dy_meters)
    
    final_gps = direction_y['lat2'], direction_x['lon2']

    print(f'\n---------------Method II-------------')
    print(f"Final GPS: {final_gps}")
    print(f'azi1_x: {lx.azi1}')
    print(f'azi1_y: {ly.azi1}')
    print(f'pixel-scale-x: {pixel_scale_x}')
    print(f'pixel-scale-y: {pixel_scale_y}')
    print(f'dx in meters: {dx_meters:.2f}m')
    print(f'dy in meters: {dy_meters:.2f}m')

    return final_gps

def find_tile_using_pixel_xy(first_tile, pixel_coordxy, zoom=18):
    """
        Given a pixel coordinate, calculate the target tile and (dx, dy) for tile-level coordinate.   
        dx, dy are tile level coordinates, which are in the range 0-256 representing, x, y offset `inside` of the given tile.

    Args:
        first_tile (tuple(int,int)): The first tile coordinate in our map image (0,0)(the upper right most corner of our image)
        pixel_coord (tuple(in,int)): The pixel coordinate for which we want the gps (with respect to initial tile location)
        zoom (int): the zoom level at which the tiles are specified (Default = 18)
    """

    row_init, col_init = first_tile
    x,y = pixel_coordxy

    # Calculate the steps needed to get from initial tile
    # to get to the target tile where our pixels reside in
    # since tiles are 256x256, we simply divide by 256 and this
    # will give us the steps needed for rows and columns 
    # example: x = 1022, y = 769 
    step_x = int(x/256)
    step_y = int(y/256)
    
    # Since we need to know where these global pixel coordinate reside exactly 
    # in our 256by256 tile internal coordinate. 
    # we do this by calculating the offset by taking the remainder from the division by 256 
    dx = abs(x % 256)
    dy = abs(y % 256) 

    selected_tile_row = row_init + step_y
    selected_tile_col = col_init + step_x
    selected_tile = (selected_tile_row, selected_tile_col)

    print(f'Pixel coordinate query: {pixel_coordxy}')
    print(f"Relative offsets inside tile: dx:{int(dx)} dy:{int(dy)}")    
    print(f'First tile:{first_tile} xy-step:({step_x},{step_y})')
    print(f'Selected tile:{selected_tile}')

    return selected_tile, (dx, dy)

def find_gps_using_pixel_xy(first_tile, pixel_coordxy, zoom=18):

    selected_tile, (dx, dy) = find_tile_using_pixel_xy(first_tile=first_tile, 
                                                       pixel_coordxy=pixel_coordxy,
                                                       zoom=zoom)
    # calculate the gps coordinate
    selected_tile_gps = calculate_latitue_longitude_from_tile(selected_tile[0], selected_tile[1], zoom)

    # first method - using approximate scale for x,y
    final_gps1 = calculate_gps_from_pixel_offsets(target_tile=selected_tile, 
                                                  offset_x=dx, 
                                                  offset_y=dy,
                                                  use_offsets_for_bearing=False,
                                                  zoom = zoom)
    
    # second method(the most accurate) - using separate pixel scales for x,y
    final_gps2 = calculate_gps_from_pixel_offsets_accurate(target_tile=selected_tile, 
                                                           offset_x=dx, 
                                                           offset_y=dy)
    
    print(f'Selected GPS:      {selected_tile_gps}')
    print(f'Final GPS-I:       {final_gps1}')
    print(f'Final GPS-II(acc): {final_gps2}')
    print(f'Distance from selected tile(GPS-I):  {get_distance(selected_tile_gps, final_gps1):.2f}m.')    
    print(f'Distance from selected tile(GPS-II): {get_distance(selected_tile_gps, final_gps2):.2f}m.')    
    print(f'-----------------------------------')

    return final_gps1, final_gps2

test:

init_gps = (38.213367074438175, -98.20816040039062)  
# when caclulating the tile coordinate using gps, make sure you account 
# for the displacement. see my note for this function above
tile = calculate_tile_coordinate(init_gps[0], init_gps[1], 18) # (100919, 59559)
tile_gps = calculate_latitue_longitude_from_tile(tile[0], tile[1], 18)
# they should match and indeed they do!
print(f'tile {tile} gps:{init_gps}')
print(f'tile {tile} gps:{tile_gps}')
# now lets test
# get the gps location at pixel coordinate (50,100)
# remember this pixel coordinate is relative to our first_tile
# we assumed that pixel(0,0) of our image is where our x0,y0
# of our first tile is.
find_gps_using_pixel_xy(tile, (50,100), 18)

output:

tile (100919, 59559) gps:(38.213367074438175, -98.20816040039062)
tile (100919, 59559) gps:(38.213367074438175, -98.20816040039062)
Pixel coordinate query: (50, 100)
Relative offsets inside tile: dx:50 dy:100
First tile:(100919, 59559) xy-step:(0,0)
Selected tile:(100919, 59559)

---------------Method I-------------
distance offset in x axis: 23.46m.
distance offset in y axis: 46.92m.
GPS based on dsx: (38.21336707, -98.20789252) az1: 90.00°
GPS based on dsy: (38.21294437, -98.20816040) az2: 180.00°
GPS based on dsx-dsy: (38.21294437, -98.20789252).
Tile GPS: (38.213367074438175, -98.20816040039062)
Final GPS:(38.212944374143305, -98.20789252325251)
Distance from tile gps to calculared gps location: 52.46m.

---------------Method II-------------
Final GPS: (38.21294558225481, -98.20789217948914)
azi1_x: 89.9995752467753
azi1_y: 180.0
pixel-scale-x: 0.46980161807460397
pixel-scale-y: 0.46785849030688803
dx in meters: 23.49m
dy in meters: 46.79m
--------------------------------------
Target Tile GPS:   (38.213367074438175, -98.20816040039062)
Final GPS-I:       (38.212944374143305, -98.20789252325251)
Final GPS-II(acc): (38.21294558225481, -98.20789217948914)
Distance from selected tile(GPS-I):  52.46m.
Distance from selected tile(GPS-II): 52.35m.