After couple of days of searching and testing I finally managed to get this to work.
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)
Also, there is this gist, where it calculates the direction of a line, in the same manner, suggested here.
# 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
Apart from this, I found out that calculating the bearing in my case can be done in a relatively simple way.
Since The azimuth is the heading measured clockwise from north, and our image is already aligned toward north, the point to the right
will have 90 degree (clockwise) angle and the point below will have
180 degree angle.
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)
Before I continue, I want to talk about the pixel scale calculation.
without a precise pixel scale calculation everything is in vein.
Searching for how to calculate the pixel scale ends up in lots of results that dont work! like the one explained here :
# 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);
The thing is, if you use the snippet above, you'll get huge deviations based on where in the world you use it to calculate the pixel scale.
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
Now, somehow in order to get this, and calculate the next latitude/longitude, I need to have the latitude of that location! it seems like a circle now!
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
Now that we have the tile coordinate, we can check at which tile our desired pixel resides, then we get that tile's gps coordinate, now we can use it's latitude to calculate the pixel scale, and then calculate the offset inside that tile, to get the exact gps location of specified pixel coordinate.
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).
This is my quick and dirty code that I have been playing wit,but before that, for better visualization this might help:
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.
xandydefined and in what reference frame are they given? You need to know.