2

I want to find the area in square metre for a polygon created from latitude and longitude.

import shapely.ops as ops
from shapely.geometry import Polygon
from pyproj import transform, Proj
from functools import partial

value = [
    [83.3203125, 58.26328705248601],
    [98.7890625, 58.81374171570782],
    [105.1171875, 60.930432202923335],
    [104.0625, 65.6582745198266],
    [97.734375, 67.47492238478702],
    [87.890625, 67.06743335108298],
    [79.8046875, 65.36683689226321],
    [79.1015625, 62.431074232920906],
    [83.3203125, 58.26328705248601],
] # from geojson

polygon = Polygon(value)

print(polygon.area)  # 191.56938242734225

The value 191.56938242734225 is not what I want, so I did some search online and found out I had to transform and use pyproj.

area = ops.transform(
    partial(transform, Proj("EPSG:4326"), Proj(proj="aea", lat_1=polygon.bounds[1], lon_1=polygon.bounds[3])), polygon
)

print(area.area)  # nan

But I get nan, what am I doing wrong here? As per a comment on the link above, for some of the data where I do get this working (for a different latitude longitude list) it does not match with the area shown in the image.

How can I get the area as seen in this image (from geojson.io)? enter image description here

Edit (after trying the solution)

from pyproj import Geod
from shapely.geometry import Polygon
from shapely.ops import orient

value = [
    (49.16848657161892, -122.4927565020954),
    (49.17145542056141, -122.44026815784645),
    (49.168944059578955, -122.42223809616848),
    (49.18182618668267, -122.42070673842014),
    (49.18901155779426, -122.40938116907655),
    (49.19333350498177, -122.41092556489616),
    (49.19754283081477, -122.4063781772052),
    (49.20350953892491, -122.38609910402336),
    (49.21439468669444, -122.36070237276812),
    (49.222360646955416, -122.35469638902534),
    (49.22411415726248, -122.35687285265936),
    (49.22498357252119, -122.35828854882728),
    (49.226469956724024, -122.35790244987238),
    (49.22832086322251, -122.35635805405282),
    (49.22936323360339, -122.35609949522244),
    (49.22913888972668, -122.35764389104204),
    (49.22773671741685, -122.35871638813896),
    (49.226739672517496, -122.36025878439766),
    (49.22010672437924, -122.38006780577793),
    (49.21996732086788, -122.38057257126464),
    (49.21999536911368, -122.39299208764706),
    (49.21992686126862, -122.4036494466543),
    (49.21891709518424, -122.40613785636424),
    (49.21888904632648, -122.4073390531128),
    (49.219562214519144, -122.40965564684215),
    (49.2199548917306, -122.41051364451972),
    (49.22074023679359, -122.4168628273335),
    (49.22121704735026, -122.421925013631),
    (49.22035066674655, -122.42482694635817),
    (49.21995799267972, -122.43289212452709),
    (49.22046285876402, -122.5120853101642),
    (49.219966402618745, -122.55650262699096),
    (49.255967240082846, -122.55753222420404),
    (49.256303572742965, -122.59614211969338),
    (49.27154824650424, -122.59648531876442),
    (49.3008416630168, -122.59480957171624),
    (49.300169602582415, -122.6655085803457),
    (49.29456874254022, -122.66036059428048),
    (49.28941538916121, -122.6675677747718),
    (49.288389695754155, -122.68427714208296),
    (49.27998615908695, -122.6956027114265),
    (49.25981182761499, -122.7177390515071),
    (49.25420638251272, -122.72099944268174),
    (49.241760018963866, -122.73953219251663),
    (49.227520721395656, -122.77317064879188),
    (49.22633007199472, -122.7634995116452),
    (49.21729887187899, -122.75397573742444),
    (49.20512626699452, -122.70094375077176),
    (49.19625956267232, -122.6707422325223),
    (49.19968297235687, -122.65383967827474),
    (49.210563761916895, -122.62311572400029),
    (49.21028322521899, -122.60552677161068),
    (49.19894821132008, -122.58853841759532),
    (49.1876984184996, -122.58114408968196),
    (49.16741121321411, -122.52297464048324),
]

polygon = Polygon(value)
geod = Geod(ellps="WGS84")
poly_area, poly_perimeter = geod.geometry_area_perimeter(orient(polygon))

print(poly_area)

I still get nan for this.

EDIT : I had my values swapped, it must be "longitude", "latitude". Thanks to snowman2 at github

python_user
  • 5,375
  • 2
  • 13
  • 32
  • Why you do not like 191? In any case, you do not tell us which datum you are using ("lat/long" coordinates doesn't tell us much, and so libraries will not understand it). On the output, you should tell it you want meter (you you need to transform to a "meter coordinate system" (usually noted with x,y). On every step you should write a comment about "datum" and "units". This will help a lot to find problems and missing conversions. – Giacomo Catenazzi Jun 25 '21 at 07:58
  • @GiacomoCatenazzi thank you for the response, I do not want `191` as my aim here is to get the area in m^2, and as per the code in the first link in my post, the answer is supposed to be in m^2, can you explain what you mean "which datum you are using", is there a way I can tell shapely or pyproj that I am using latitude and longitude? – python_user Jun 25 '21 at 09:28
  • Latitude and longitude are just a way to call two angles (like using x,y), but it doesn't really define precise point to earth. You need to specify how you measure the angles, etc. (this is the datum). Different countries (and different times) uses different origins/offset (and also angles), but also the form of the Earth. All geographic program needs such datum (you gave one on pyproj: the WGS84 (aka the "GPS coordinates"). Datum is one of the most confusing things working with projections (and we all make a lot of errors): so write it explicit (and double check). – Giacomo Catenazzi Jun 25 '21 at 13:21
  • @GiacomoCatenazzi ok this makes sense to me, I thought this would be a task that would not require much domain knowledge, but I was wrong – python_user Jun 25 '21 at 14:22

1 Answers1

3

You could use pyproj to get the geodesic area: https://pyproj4.github.io/pyproj/stable/examples.html#geodesic-area

from pyproj import Geod
from shapely.geometry import Polygon


value = [
    [83.3203125, 58.26328705248601],
    [98.7890625, 58.81374171570782],
    [105.1171875, 60.930432202923335],
    [104.0625, 65.6582745198266],
    [97.734375, 67.47492238478702],
    [87.890625, 67.06743335108298],
    [79.8046875, 65.36683689226321],
    [79.1015625, 62.431074232920906],
    [83.3203125, 58.26328705248601],
]

polygon = Polygon(value)
geod = Geod(ellps="WGS84")
poly_area, poly_perimeter = geod.geometry_area_perimeter(polygon)

print(poly_area)  # 1073367471345.5327

Also see: https://pyproj4.github.io/pyproj/stable/api/geod.html#pyproj.Geod.geometry_area_perimeter

You may need shapely.ops.orient

snowman2
  • 646
  • 4
  • 11
  • so what is the value of `geod`? the code given here does not have that, from the link `geod = Geod('+a=6378137 +f=0.0033528106647475126')` can you explain or point me what `'+a=6378137 +f=0.0033528106647475126'` is and should I change it? – python_user Jun 25 '21 at 14:17
  • I get `nan` for a different data (assume, data_1) and I get results in negative for a different data (assume, data_2), I will add data_1 and data_2 when I can access the code, do you mind explaining that? I used `geod` value `'+a=6378137 +f=0.0033528106647475126'` – python_user Jun 25 '21 at 14:28
  • Whoops, forgot geod. Added it in. – snowman2 Jun 25 '21 at 14:33
  • appreciate the quick edit, can you try with https://pastebin.com/wuc5ZB2L, that has the `data_1` and `data_2`, which give me a `nan` and negative area – python_user Jun 25 '21 at 14:38
  • Did you use `shapely.ops.orient` for the negative area? – snowman2 Jun 25 '21 at 15:11
  • no I have not, can you show an example if you dont mind? – python_user Jun 25 '21 at 15:58
  • You just do it on the polygon before passing it into the function. – snowman2 Jun 25 '21 at 16:22
  • negative area is no more after `orient`, but I still get `nan`, any workaround for that? I have edited the code I tried – python_user Jun 28 '21 at 02:23
  • Not sure about the nan. Might be worth raising an issue: https://github.com/pyproj4/pyproj – snowman2 Jun 28 '21 at 03:02