0

I would like to create a polygon assuming it is drawn in a geographic environment using latitude, longitude, bearing, length and width. The code is below:

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

The variables are:

v1_coord_5= [41.39129840757358, 2.1658667401858738]
v1_length = 5 #meters
v1_width = 1.8 #meters
v1_bearing=148.26196928529635

points = get_rotated_points(tuple(v1_coord_5), v1_bearing,
                               v1_width, v1_length)

The points results are the following:

POLYGON ((41.39127237180109 2.1658768035733806, 41.39128615544415 2.1658881248848436, 41.391324443338675 2.165856676778987, 41.39131065969457 2.1658453554632664, 41.39127237180109 2.1658768035733806))

So far, I've tried the following: 1)

!pip install area
from area import area
obj = {'type':'Polygon','coordinates':[[[41.39127237180109, 2.1658768035733806],[41.39128615544415, 2.1658881248848436],[41.391324443338675, 2.165856676778987],[41.39131065969457, 2.1658453554632664],[41.39127237180109, 2.1658768035733806]]]}
area(obj)

Result: 10.735455561073895

2)

!pip install pyproj

co = {"type": "Polygon", "coordinates": [
    [(41.39127237180109, 2.1658768035733806),
     (41.39128615544415, 2.1658881248848436),
     (41.391324443338675, 2.165856676778987),
     (41.39131065969457, 2.1658453554632664)]]}
lon, lat = zip(*co['coordinates'][0])
from pyproj import Proj
pa = Proj("+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55")
x, y = pa(lon, lat)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
from shapely.geometry import shape
shape(cop).area  # 268952044107.43506

Result: 10.663792266425144

3)

from shapely.geometry import Polygon

points = [(41.39127237180109, 2.1658768035733806), (41.39128615544415, 2.1658881248848436), (41.391324443338675, 2.165856676778987), (41.39131065969457, 2.1658453554632664), (41.39127237180109, 2.1658768035733806)]
polygon = Polygon(points)
area_sdeg = polygon.area

Result: 8.669387666688732e-10

If I use a regular formula, the result is 9 m2. Unfortunately, I couldn't get that result. Any clue what I'm missing? thanks

GoT GOt
  • 49
  • 1
  • 9
  • Looks like duplicate of: https://stackoverflow.com/questions/4681737/how-to-calculate-the-area-of-a-polygon-on-the-earths-surface-using-python – Gameplay Dec 02 '22 at 07:24
  • Which "regular formula"? In any case there is no true answer. The area depends on the shape of the geoid you use (and there are many possible choices), coordinate system you use (yeah: there are different way you can measure latitude, and 41 degree is in the "bad" place where it matter more), etc. And in theory there are two areas (usually we want the small one (the sum of the two is surface of Earth), in any case often you should define the polygon with the right direction. Without defining the base geometry no value are meangful – Giacomo Catenazzi Dec 02 '22 at 07:50

0 Answers0