What is the cause for the degradation of environment?
Capitalism, corruption, consuming society? - OVERPOPULATION!
Please, save the Planet - kill yourself...

Friday, September 21, 2012

Utility to Calculate Polygon Width in Any Given Direction

I have created an "Azimuth-Width" utility recently. It is designed to compute (for now) maximum or minimum width of the given polygon(s) in the given direction. Corresponding QGIS plugin coming soon (or not too soon), so this utility may help these who are in need right now.


Prereqirements:

  1. You have to have QGIS installed on your machine.
  2. QGIS have to be in PYTHONPATH. See: http://qgis.org/pyqgis-cookbook/intro.html
  3. If you are non-Windows user ensure that command "QgsApplication.setPrefixPath( qgis_prefix, True)" in the utility file have  a valid path to QGIS installation. If not - modify "qgis_prefix" to set correct value. It is "/usr" now - must work in most of cases.

General usage:

Copy width.py to the directory with a shh-file containing polygons.

Using console navigate to the directory. In console type:
:~> python ./width.py [file to analyse] [field to store values] [azimuth (decimal degrees)] [mode ('min' or 'max')] [mode-2 ('abs', or 'rel')] [algorithm ('byStep' or 'byVertex' or 'Mix')] [step (real numver, CRS units; for 'byStep' and 'Mix' only)]

Example:~> python ./width.py poly.shp width 285.9 max abs byStep 1.3

Where:
  • ./width.py - name of this utility file.
  • [file to analyse] - a shp-file with polygons to analyse.
  • [field to store values] - if it does not exist it will be created.
  • [azimuth] - a direction for width calculation. Accepts decimal degrees from 0 to 360.
  • [mode] - type of width to calculate. Currently 'min' (returns minimum width value in given direction) and 'max' (returns maximum value in given direction) modes are available. Mode 'min' used with 'byVertex' algorithm will return minimum polygon width different then 0.0. This will be different (greatly most time) to 'min' mode and 'byStep' or 'Mix' algorithm.
  • [mode-2] - if polygon is not convex it may have several segments in given direction. If you want to take sum of the segments of the result - use 'abs', if you want only the longest (shortest) segment - use 'rel'. Fixed: Note that currently 'rel' sometimes will provide incorrect results for non-convex polygons and for convex polygons with redundant vertices on its edges. This issue will be solved when the behaviour of the intersection() command will be changed or when I will implement a workaround for it.
  • [algorithm] - algorithm that will be used: 'byStep', 'byVertex' or 'Mix'. "byStep" algorithm (improved version of the general idea that described here) will take provided value (shp-file CRS units) and swipe the polygons line with the given step along the bounding box. Precision depends on the step. Speed and precision depends on step - lower step mean more precise result but computation will take more time. "byVertex" (general idea described here) will intersect polygon with the line only in vertexes of the given polygon. For convex polygons this algorithm will be faster and more precise using 'max' mode then "byStep" algorithm. But this algorithm is not suitable for 'min' mode. "Mix" algorithm will use both "byVertex" and "byStep" algorithm so in some cases it will be most precise but will consume even more time than 'byStep'.
  • [step] - step for "byStep" and "Mix" algorithm. Takes decimal values in shp-file CRS units. Lower step means more precision but more computation time.

Some advises:

  1. Use Equal Area projections.
  2. Make sure your polygons have correct geometry. Otherwise the calculations will be incorrect.
  3. 'byStep' algorithm should be faster when polygon have enormous number of vertices.
  4. If you have a large variety in polygons area e.g. a continent and an island to save computation time you may define big step (like 1000 m or so) to swipe through continent faster. It will save computation time and a width for the small island will be calculated even if step is grater then any side of island's' Bounding Box: step for the island in this case will be 1/100 of the shortest side of it's Bounding Box.

1 comment :