Search code examples
dictionarycoordinateslatitude-longitudefoliumshapely

How to calculate lat and lon of a rectangle from a center point


I would like to draw a rectangle based on a center point lat and lon assuming a given length and width, let's say 4.5m and 1.5m, respectively. I guess, we need the bearing too. I've made a simulation by drawing a rectangle on Google Earth, getting the positions and putting them on my code. However, I need something automatic. My question is how can I link the Cartesian coordinates to those four points (rectangle) in meters.

Here is my code:

import geopandas as gpd
from shapely.geometry import Polygon

lat_point_list = [41.404928,  41.404936,  41.404951,  41.404943]
lon_point_list = [2.177339,   2.177331,   2.177353,   2.177365]

polygon_geom = Polygon(zip(lon_point_list, lat_point_list))
import folium
m = folium.Map([41.4049364, 2.1773560], zoom_start=20)
folium.GeoJson(polygon_geom).add_to(m)
folium.LatLngPopup().add_to(m)
m

I would like this:

enter image description here

Update:

I know this is basic trigonometry. If I split the rectsngle into triangles, we can find the different points. I know it is basic for simple exercises, however, I don't know of it changes when using Cartesian coordinates. Then, my goal is to get the points A, B, C and D, knowing the center of the rectangle in latitude and longitude, length and width.

enter image description here


Solution

  • Get the rectangular (NE, SW) bounds of your point and use that as bounds to folium.Rectangle.

    Example, using your data. 4.5m and 1.5m are a bit small to see the rectangle:

    import geopy
    import geopy.distance
    import math
    import folium
    
    def get_rectangle_bounds(coordinates, width, length):
        start = geopy.Point(coordinates)
        hypotenuse = math.hypot(width/1000, length/1000)
    
        # Edit used wrong formula to convert radians to degrees, use math builtin function
        northeast_angle = 0 - math.degrees(math.atan(width/length)) 
        southwest_angle = 180 - math.degrees(math.atan(width/length)) 
    
        d = geopy.distance.distance(kilometers=hypotenuse/2)
        northeast = d.destination(point=start, bearing=northeast_angle)
        southwest = d.destination(point=start, bearing=southwest_angle)
        bounds = []
        for point in [northeast, southwest]:
            coords = (point.latitude, point.longitude)
            bounds.append(coords)
    
        return bounds
    
    # To get a rotated rectangle at a bearing, you need to get the points of the the recatangle at that bearing
    def get_rotated_points(coordinates, bearing, width, length):
        start = geopy.Point(coordinates)
        width  = width/1000
        length = length/1000
        rectlength = geopy.distance.distance(kilometers=length)
        rectwidth = geopy.distance.distance(kilometers=width)
        halfwidth = geopy.distance.distance(kilometers=width/2)
        halflength = geopy.distance.distance(kilometers=length/2)
    
        pointAB = halflength.destination(point=start, bearing=bearing)
        pointA = halfwidth.destination(point=pointAB, bearing=0-bearing)
        pointB = rectwidth.destination(point=pointA, bearing=180-bearing)
        pointC = rectlength.destination(point=pointB, bearing=bearing-180)
        pointD = rectwidth.destination(point=pointC, bearing=0-bearing)
    
        points = []
        for point in [pointA, pointB, pointC, pointD]:
            coords = (point.latitude, point.longitude)
            points.append(coords)
    
        return points
    
    start_coords = [41.4049364, 2.1773560]
    length = 4.50 #in meters
    width = 1.50
    bearing = 45 #degrees
    
    m = folium.Map(start_coords, zoom_start=20)
    bounds = get_rectangle_bounds(tuple(start_coords),width, length )
    points = get_rotated_points(tuple(start_coords), bearing, width, length)
    
    folium.Rectangle(bounds=bounds,
                    fill=True,
                    color='orange',
                    tooltip='this is Rectangle'
                   ).add_to(m)
    
    # To draw a rotated rectangle, use folium.Polygon
    folium.Polygon(points).add_to(m)