Skip to content

Using corners of the map to calculate the area of whole world does not return correct result #1074

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
masterAllen opened this issue May 9, 2022 · 5 comments
Labels

Comments

@masterAllen
Copy link

Code Sample, a copy-pastable example if possible

from pyproj import Geod
from shapely.ops import transform, orient
from shapely.geometry import Point, Polygon

data = [
    (89.9999, 179.9999),
    (89.9999, -179.9999),
    (-89.9999, -179.9999),
    (-89.9999, 179.9999),
]
world_poly = Polygon(orient([(y, x) for (x, y) in data]))

geod = Geod(ellps="WGS84") 
area = abs(geod.geometry_area_perimeter(world_poly)[0])
print(area)

Problem description

I want to calculate the area of whole world, so i use polygon with four points which are corners of the map. But the result is 283369789.8556872 which is very small.

Expected Output

The area of whole world

Environment Information

  • Output from: pyproj -v
pyproj info:
    pyproj: 3.0.1
      PROJ: 8.0.1
  data dir: /usr/share/proj
user_data_dir: /home/allen/.local/share/proj

System:
    python: 3.10.4 (main, Mar 23 2022, 23:05:40) [GCC 11.2.0]
executable: /usr/bin/python
   machine: Linux-5.17.4-arch1-1-x86_64-with-glibc2.35

Python deps:
       pip: 21.0
setuptools: 59.5.0
    Cython: 0.29.28

Installation method

  • pip
@masterAllen masterAllen added the bug label May 9, 2022
@snowman2
Copy link
Member

@cffk, apologies for bothering you again. However, I think you would be able to give a much better response to this question.

@cffk
Copy link
Contributor

cffk commented May 12, 2022

Plot your 4 plots on a globe instead of a map (and replace the 89.9999 by 90) and you will see that they outline a narrow lune of width 0.0002 deg. If you multiply your result by 360/0.0002 you get the full area of the earth. You got a negative result because your lune was traversed clockwise.

Better would be to bite off half the globe, e.g., with 3 points on the equator

  data = [ (0, 0), (0, 120), (0, -120) ]

and multiply this result by 2.

geometry_area_perimeter won't return more than half the area of the earth. A larger polygon will be treated as enclosing the rest of the earth in the opposite direction. This is the issue that you ran into.

@cffk
Copy link
Contributor

cffk commented May 15, 2022

@masterAllen Are you set now? Unless you still have questions, this issue can be closed as resolved.

@snowman2 snowman2 added question and removed bug labels May 15, 2022
@snowman2
Copy link
Member

Thanks @cffk 👍. @masterAllen feel free to comment if you have further questions.

@masterAllen
Copy link
Author

So sorry for late reply. I have no questions now, thank you for your help :D

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants