[GIS] Calculate vertices’ coordinates of a rectangle, given midline and broadth of rectangle

coordinate systemcoordinatesextentsgeographyvertices

I have latitude and longitude of two points making the midline of a rectangle. Now I want to determine coordinates of 4 vertices of that rectangle. I know how to calculate on plane geometry by using line equation. But on geographic geometry, I'm new. I already read how to calculate one unit of latitude and longitude decimal degree but I'm still confused how to add them to a given point to calculate wanted points. My problem is like this:

Given:

Broadth of rectangle
Coordinates of two points making the midline, distance between these two points are equal to the width

Result: Coordinates of 4 vertices of rectangle

For example:
The given midline has two points:

A (10.767008,106.665884)
B (10.767715,106.667151)

Broadth is 4-metre long

All I could do now is to calculate the width of rectangle from those two points. How can I use width and height to calculate 4 vertices of this rectangle?

All of my calculated cases are on small scale (part of one street of a city) so I think altitude can be excluded

Best Answer

Being that you are dealing with such a small spatial area, we can get the four coordinate pairs using a couple of Euclidean equations. Unfortunately, GIS.SE doesn't support LaTeX equations, so I am going to assume you know basic geometry equations like slope, point-slope form, and the euclidean distance formula.

To better illustrate what is being done, refer to this image:

The goal is to find the coordinates of r, s, t, and u.

First, let's state the givens.

A (lat: 10.767008, lng: 106.665884) 
B (lat: 10.767715, lng: 106.667151)
r-s, t-u = 4 metres or roughly .000036036 decimal degrees at this latitude

(if you need a more accurate way of determining this, there are plenty of posts on the subject)

Second, calculate the perpendicular slope of line A-B, we call it M.

Similar to slope equation, however with perpendicular slope your numerator is:

x_1 - x_2

And your denomimator is:

y_2 - y_1

This means:

M = (106.665884 - 106.667151)/(10.767715 - 10.767008)
M = -1.7920792

Third, we write out our equations for our segments, solving for the y coordinates:

y_r = M(x_r - 106.667151) + 10.767715
y_s = M(x_s - 106.667151) + 10.767715
y_t = M(x_t - 106.665884) + 10.767008
y_u = M(x_u - 106.665884) + 10.767008

Fourth, accounting for slope being positive or negative, we then write out the equations solving for the x coordinate of our unknown points. If the slope is negative, the x coordinate of s and u will be the result of the negative difference. If positive, the x coordinate of s and u will be the positive sum. Since our slope is negative, that means:

x_s = 106.667151 - (0.000018018 / sqrt(1 + (-1.7920792)^2))
x_s = 106.667139

x_r = 106.667151 + (0.000018018 / sqrt(1 + (-1.7920792)^2))
x_r = 106.667163

x_t = 106.665884 + (0.000018018 / sqrt(1 + (-1.7920792)^2))
x_t = 106.665896

x_u = 106.665884 - (0.000018018 / sqrt(1 + (-1.7920792)^2))
x_u = 106.665872

Fifth, now that we have our x coordinates, we can plug the values back into our segment equations and get our y coordinates.

y_r = -1.7920792 * (106.667163 - 106.667151) + 10.767715
y_r = 10.7676935

y_s = -1.7920792 * (106.667139 - 106.667151) + 10.767715
y_s = 10.7677365

y_t = -1.7920792 * (106.665896 - 106.665884) + 10.767008
y_t = 10.7669865

y_u = -1.7920792 * (106.665872 - 106.665884) + 10.767008
y_u = 10.7670295

So the coordinates to create your polygon would be:

r (lat: 10.7676935, lng: 106.667163)
s (lat: 10.7677365, lng: 106.667139)
t (lat: 10.7669865, lng: 106.665896)
u (lat: 10.7670295, lng: 106.665872)

Understand, they may be slightly off due to the assumption of metre to decimal degree length at this latitude, and values were rounded off.

Related Solutions

[GIS] Getting all vertex lat, long coordinates every 1 meter between two known points using PHP

This is quick and dirty. Most of this was adapted from the following website:

http://www.movable-type.co.uk/scripts/latlong.html

I haven't tested this too much but it seemed to work after initial testing. It will return a list of lat, long coordinate pairs along a line at a specified interval. I wrote it in python since I don't know php very well.

from math import cos, sin, atan2, sqrt, radians, degrees, asin, modf

def getPathLength(lat1,lng1,lat2,lng2):
    '''calculates the distance between two lat, long coordinate pairs'''
    R = 6371000 # radius of earth in m
    lat1rads = radians(lat1)
    lat2rads = radians(lat2)
    deltaLat = radians((lat2-lat1))
    deltaLng = radians((lng2-lng1))
    a = sin(deltaLat/2) * sin(deltaLat/2) + cos(lat1rads) * 
        cos(lat2rads) * sin(deltaLng/2) * sin(deltaLng/2)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = R * c

    return d

def getDestinationLatLong(lat,lng,azimuth,distance):
    '''returns the lat an long of destination point 
    given the start lat, long, aziuth, and distance'''
    R = 6378.1 #Radius of the Earth in km
    brng = radians(azimuth) #Bearing is degrees converted to radians.
    d = distance/1000 #Distance m converted to km

    lat1 = radians(lat) #Current dd lat point converted to radians
    lon1 = radians(lng) #Current dd long point converted to radians

    lat2 = asin( sin(lat1) * cos(d/R) + cos(lat1)* sin(d/R)* cos(brng))

    lon2 = lon1 + atan2(sin(brng) * sin(d/R)* cos(lat1), 
           cos(d/R)- sin(lat1)* sin(lat2))

    #convert back to degrees
    lat2 = degrees(lat2)
    lon2 = degrees(lon2)

    return[lat2, lon2]

def main(interval,azimuth,lat1,lng1,lat2,lng2):
    '''returns every coordinate pair inbetween two coordinate 
    pairs given the desired interval'''
    coords = []
    d = getPathLength(lat1,lng1,lat2,lng2)
    remainder, dist = modf((d / interval))
    counter = 1.0
    coords.append([lat1,lng1])
    for distance in xrange(1,int(dist)):
        c = getDestinationLatLong(lat1,lng1,azimuth,counter)
        counter += 1.0
        coords.append(c)
    counter +=1
    coords.append([lat2,lng2])
    return coords

if __name__ == "__main__":
    #point interval in meters
    interval = 1
    #direction of line in degrees
    azimuth = 279.91
    #start point
    lat1 = 41.97809
    lng1 = -91.65573
    #end point
    lat2 = 41.97834
    lng2 = -91.65767
    coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
    print coords

I hope this helps. sorry if not. :)

Ok, I added the ability for it to calculate the azimuth for you. I tested with the same location you were using and it turned out perfectly.http://www.darrinward.com/lat-long/?id=705140

import math

def getPathLength(lat1,lng1,lat2,lng2):
    '''calculates the distance between two lat, long coordinate pairs'''
    R = 6371000 # radius of earth in m
    lat1rads = math.radians(lat1)
    lat2rads = math.radians(lat2)
    deltaLat = math.radians((lat2-lat1))
    deltaLng = math.radians((lng2-lng1))
    a = math.sin(deltaLat/2) * math.sin(deltaLat/2) + math.cos(lat1rads) * math.cos(lat2rads) * math.sin(deltaLng/2) * math.sin(deltaLng/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = R * c
    return d

def getDestinationLatLong(lat,lng,azimuth,distance):
    '''returns the lat an long of destination point 
    given the start lat, long, aziuth, and distance'''
    R = 6378.1 #Radius of the Earth in km
    brng = math.radians(azimuth) #Bearing is degrees converted to radians.
    d = distance/1000 #Distance m converted to km
    lat1 = math.radians(lat) #Current dd lat point converted to radians
    lon1 = math.radians(lng) #Current dd long point converted to radians
    lat2 = math.asin(math.sin(lat1) * math.cos(d/R) + math.cos(lat1)* math.sin(d/R)* math.cos(brng))
    lon2 = lon1 + math.atan2(math.sin(brng) * math.sin(d/R)* math.cos(lat1), math.cos(d/R)- math.sin(lat1)* math.sin(lat2))
    #convert back to degrees
    lat2 = math.degrees(lat2)
    lon2 = math.degrees(lon2)
    return[lat2, lon2]

def calculateBearing(lat1,lng1,lat2,lng2):
    '''calculates the azimuth in degrees from start point to end point'''
    startLat = math.radians(lat1)
    startLong = math.radians(lng1)
    endLat = math.radians(lat2)
    endLong = math.radians(lng2)
    dLong = endLong - startLong
    dPhi = math.log(math.tan(endLat/2.0+math.pi/4.0)/math.tan(startLat/2.0+math.pi/4.0))
    if abs(dLong) > math.pi:
         if dLong > 0.0:
             dLong = -(2.0 * math.pi - dLong)
         else:
             dLong = (2.0 * math.pi + dLong)
    bearing = (math.degrees(math.atan2(dLong, dPhi)) + 360.0) % 360.0;
    return bearing

def main(interval,azimuth,lat1,lng1,lat2,lng2):
    '''returns every coordinate pair inbetween two coordinate 
    pairs given the desired interval'''

    d = getPathLength(lat1,lng1,lat2,lng2)
    remainder, dist = math.modf((d / interval))
    counter = float(interval)
    coords = []
    coords.append([lat1,lng1])
    for distance in xrange(0,int(dist)):
        coord = getDestinationLatLong(lat1,lng1,azimuth,counter)
        counter = counter + float(interval)
        coords.append(coord)
    coords.append([lat2,lng2])
    return coords

if __name__ == "__main__":
    #point interval in meters
    interval = 1.0
    #direction of line in degrees
    #start point
    lat1 = 43.97076
    lng1 = 12.72543
    #end point
    lat2 = 43.969730
    lng2 = 12.728294
    azimuth = calculateBearing(lat1,lng1,lat2,lng2)
    print azimuth
    coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
    print coords

I think I have fixed the bugs from the original and it should work (including change the interval). Let me know how it turns out.

[GIS] Longitude difference to Meters

Just in case you want the distance along the parallel of latitude (and not the shortest distance) between two points at a given latitude,

distance = r * (longitude 2 - longitude 1)

with longitude in radians and where r, the radius of the parallel of latitude, is

r = R * cosine latitude

where R is Earth radius.

Best Answer

Related Solutions

[GIS] Getting all vertex lat, long coordinates every 1 meter between two known points using PHP

[GIS] Longitude difference to Meters

Related Question