Python / shapely: Calculate Polygon area in planar units (e.g. square-meters) -

i using python 3.4 , shapely 1.3.2 create polygon object out of list of long/lat coordinate pairs transform well-known-text string in order parse them. such polygon might like:

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

since shapely not handle projections , implements geometry objects in carthesian space, calling area method on polygon like:

poly.area 

gives me area of polygon in unit of square-degrees. area in planar unit square-meters, guess have transform coordinates of polygon using different projection (which one?).

i read several times pyproj library should provide way this. using pyproj, there way transform whole shapely polygon object projection , calculate area?

i other stuff polygons (not think now) , in cases, need calculate area.

so far, found example: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/

which mean splitting each polygon object outer and, if present, inner rings, grab coordinates, transform each pair of coordinates projection , rebuild polygon object, calculate area (what unit anyway?). looks solution, not practical.

any better ideas?

ok, made basemap toolkit of matplotlib library. explain how works, maybe helpful sometime.

1. download , install matplotlib library on system. http://matplotlib.org/downloads.html

for windows binaries, recommend using page: http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib beware of hint says:

requires numpy, dateutil, pytz, pyparsing, six, , optionally pillow, pycairo, tornado, wxpython, pyside, pyqt, ghostscript, miktex, ffmpeg, mencoder, avconv, or imagemagick.

hence, if not installed on system, have download , install numpy, dateutil, pytz, pyparsing , 6 matplotlib work (for windows users: of them can found on page, linux users, python package manager "pip" should trick).

2. download , install "basemap" toolkit matplotlib. either http://matplotlib.org/basemap/users/installing.html or windows binaries here: http://www.lfd.uci.edu/~gohlke/pythonlibs/#basemap

3. projection in python code:

#import necessary libraries mpl_toolkits.basemap import basemap import matplotlib.pyplot plt import numpy np  #coordinates transformed long = -112.367 lat = 41.013  #create basemap projection. 1 use you. #some examples can found @ http://matplotlib.org/basemap/users/mapsetup.html m = basemap(projection='sinu',lon_0=0,resolution='c')  #project coordinates: projected_coordinates = m(long,lat) 

output:

projected_coordinates (10587117.191355567, 14567974.064658936)

simple that. now, when using projected coordinates build polygon shapely , calculating area via shapely's area method, you'll area in unit of square-meters (according projection used). square-kilometers, divide 10^6.

edit: struggled hard not transform single coordinates, whole geometry objects polygons when given shapely objects , not via raw coordinates. meant writing lot of code to

• get coordinates of outer ring of polygon
• determine if polygon has holes , if so, process each hole seperately
• transform each pair of coordinates of outer ring , holes
• put whole thing , create polygon object projected coordinates
• and that's polygons... add more loops multipolygons , geometrycollections

then stumbled upon part of shapely documentation makes things lot easier: http://toblerity.org/shapely/manual.html#shapely.ops.transform

when projection map set, example done above:

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

then, 1 can transform shapely geometry object using projection via:

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