Calculate polygon area in flat units (e.g. square meters) in Shapely

Question

I use Python 3.4 and shapely 1.3.2 to create a Polygon object from a list of long / lat-coordinate pairs, which I convert to a well-known text string for parsing them. Such a polygon might look like this:

POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407, -116.908 43.375, -116.904 43.371)) 

Since shapely does not process any projections and implements all geometry objects in the carcinoma space, calling the region method on this polygon, as:

 poly.area 

gives me the area of ​​this polygon squared square. To get the area in a flat block, such as square meters, I assume that I would have to convert the coordinates of the polygon using a different projection (which?).

I read several times that the pyproj library should provide a way to do this. Using pyproj, is there a way to convert a whole slender Polygon object to another projector and then calculate the area?

I do some other things with my polygons (not what you think now), and only in some cases I need to calculate the area.

So far, I have found only this example: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/

which would mean splitting each Polygon object into its outer and, if any, inner rings, grab the coordinates, convert each pair of coordinates to a different projection and rebuild the Polygon object, and then calculate its area (which unit is it anyway?). This seems like a solution, but not very practical.

Any better ideas?

Answer 1

Ok, I finally did this with the base matplotlib library toolkit. I will explain how it works, maybe it will be useful to someone someday.

1. Download and install the matplotlib library on your system. http://matplotlib.org/downloads.html

For Windows binaries, I recommend using this page: http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib Beware of the hint that reads:

Requires numpy, dateutil, pytz, pyparsing, six and possibly a pillow, pycairo, tornado, wxpython, pyside, pyqt, ghostscript, miktex, ffmpeg, mencoder, avconv or imagemagick.

Therefore, if you have not already installed on your system, you need to download and install numpy, dateutil, pytz, pyparsing and six, as well as for matplotlib to work correctly (for Windows users: all of them can be found on the page, for Linux users, the package manager python "pip" should do the trick).

2. Download and install the basemap toolkit for matplotlib. Either from http://matplotlib.org/basemap/users/installing.html or for Windows binaries also from here: http://www.lfd.uci.edu/~gohlke/pythonlibs/#basemap

3. Make a projection in Python code:

 #Import necessary libraries from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt import numpy as np #Coordinates that are to be transformed long = -112.367 lat = 41.013 #Create a basemap for your projection. Which one you use is up to you. #Some examples can be found at http://matplotlib.org/basemap/users/mapsetup.html m = Basemap(projection='sinu',lon_0=0,resolution='c') #Project the coordinates: projected_coordinates = m(long,lat) 

Output:

projected_coordinates (10587117.191355567, 14567974.064658936)

Just. Now, when you use the projected coordinates to build a polygon with a beautiful pattern, and then calculate the area using the slender area method, you will get the area in square meter (in accordance with the projection you use). To get square kilometers, divide them by 10 ^ 6.

Edit: I did my best not to convert only single coordinates, but whole geometry objects, such as polygons, when they were already set as beautiful objects, and not through their original coordinates. This meant writing a lot of code in

• get the coordinates of the outer ring of the polygon
• determine if the polygon has holes, and if so, treat each hole separately
• transform each pair of coordinates of the outer ring and any holes
• combine all this and create a polygon object with projected coordinates
• and that only for polygons ... add even more loops for multipolygon and geometric collections.

Then I came across this piece of beautiful documentation that makes things a lot easier: http://toblerity.org/shapely/manual.html#shapely.ops.transform

When a projection map is installed, for example, as done above:

m = Basemap(width=1,height=1, resolution='l',projection='laea',lat_ts=50,lat_0=50,lon_0=-107.)

Then, using this projection, you can transform any slender geometry object through:

 from shapely.ops import transform projected_geometry = transform(m,geometry_object) 

Answer 2

Convert to radians and assume that the Earth is an ideal sphere of radius 6370Km:

<p> p = np.array ([[[- 116.904,43.371], [- 116.823, 43.389], [- 116.895,43.407], [- 116.908,43.375], [- 116.904,43.371]])

poly = Polygon (np.radians (p))

poly.area = 4.468737548271707e-07

poly.area * 6370 ** 2 = +18.132751662246623