Python: Contour features in a polygon shapefiles with a straightforward solution

3
votes

Dear Members list,

First of all, I apologize to post this question modified and improved from a previous post. Recently i am working on shapefiles polygon in order to compute basic contour features:

  1. Area,
  2. Perimeter,
  3. Area convex hull,
  4. Perimeter convex hull,
  5. major axis length = gives the length of major axis,
  6. minor axis length = gives the length of minor axis,

where major and minor axis length are computed following the Figure:

enter image description here

Using osgeo.gdal, ogr and shapely is possible to load and calculate all indeces but not the major and minor axis length. Reading online solution can be using

  1. scikit-image = Measure region properties
  2. OpenCV

I am looking a straightforward solution in order to make my code easy and elegant. Some blogs suggest to make an ellipse approximation to my polygon in order to retrieve the major and minor axis length. Is it the best solution?

Any references would be quite helpful. Thanks in advance 

import osgeo.gdal, ogr
from shapely.geometry import Polygon

shp = osgeo.ogr.Open('../examples/mypoly.shp')
layer = shp.GetLayer()
feature = layer.GetFeature(0)
geometry = feature.GetGeometryRef()
# get area
Area = geometry.GetArea()
pts = geometry.GetGeometryRef(0)
points = []
for p in range(pts.GetPointCount()):
   points.append((pts.GetX(p), pts.GetY(p)))
polygon = Polygon(points)
# get Perimeter
Perimeter = polygon.length
# convex Hull
ConvexHull = polygon.convex_hull
# get Perimeter convex Hull
PerimeterConvexHull = ConvexHull.length
# get Area convex Hull
AreaConvexHull = ConvexHull.area

these are the coordinate vertices of my polygon

polygon = Polygon([(560023.4495758876400000 6362057.3904932579000000),(560023.4495758876400000 6362060.3904932579000000),(560024.4495758876400000 6362063.3904932579000000),(560026.9495758876400000 6362068.3904932579000000),(560028.4495758876400000 6362069.8904932579000000),(560034.9495758876400000 6362071.8904932579000000),(560036.4495758876400000 6362071.8904932579000000),(560037.4495758876400000 6362070.3904932579000000),(560037.4495758876400000 6362064.8904932579000000),(560036.4495758876400000 6362063.3904932579000000),(560034.9495758876400000 6362061.3904932579000000),(560026.9495758876400000 6362057.8904932579000000),(560025.4495758876400000 6362057.3904932579000000),(560023.4495758876400000 6362057.3904932579000000)])

in order to test my code from here:

polygon = Polygon(points)
# get Perimeter
Perimeter = polygon.length
# convex Hull
ConvexHull = polygon.convex_hull
# get Perimeter convex Hull
PerimeterConvexHull = ConvexHull.length
# get Area convex Hull
AreaConvexHull = ConvexHull.area
1
pythonmathgeometrypolygonellipse
How exactly do you define the major and minor axis for a polygon? – ev-br
@Zhenya. Thanks to replay and for the observation. The major and minor axis are define as the same approach of openCV and MATLAB Image Processing Toolbox (regionprops) www-rohan.sdsu.edu/doc/matlab/toolbox/images/regionprops.html – Gianni Spear
the openCV link you've linked lists majoraxis_length and minoraxis_length, can you use these then? – ev-br
@Zhenya, thanks for suggestion but i wish to work with polygon (if i can). The first way, following this post gis.stackexchange.com/questions/22895/…, is calculate the minimum-area-rectangle for given points. – Gianni Spear
@Zhenya: see this new post stackoverflow.com/questions/13542855/… – Gianni Spear

1 Answers

0
votes

Hi，i have some work on shapefiles polygon likes your work. And i use the minimum_rotated_rectangle in shapely to define the majoraxis_length and minoraxis_length. and my code is it:

polygon = Polygon(points)  
ConvexHull = polygon.convex_hull  #  convex Hull
#minimum_ rotated_ Rectangle is the minimum boundary rectangle, exterior is the external point coordinates of the minimum boundary matrix, and coords is the coordinates of the points
p1 = Point(ConvexHull.minimum_rotated_rectangle.exterior.coords[0])
p2 = Point(ConvexHull.minimum_rotated_rectangle.exterior.coords[1])
p3 = Point(ConvexHull.minimum_rotated_rectangle.exterior.coords[2])
majoraxis_length  =  p1.distance(p2)
minoraxis_length = p2.distance(p3)