Pages

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.

Sunday, September 2, 2012

Raster Extent Modification Using QGIS

UPDATE: you can download script for Processing module to perform this task from here.

When you need to perform algebraic operations with two or more rasters you often have to ensure that these rasters have the same extent, otherwise the operation either  won't be allowed (in SAGA for example) or you will get output raster cropped to the overlapping part. With regular (not georeferenced) images you can simply load them in GIMP and increase or decrease canvas for amount of pixels you like. But in case of georeferenced images you will loose spatial reference information. 

In QGIS you can change extent of the rasters. Lets examine one of the worst case scenarios. There are two overlapping (one band) rasters A and B. Say, we need to add A values to B values and get the final image to have extent that will contain both images.

Overlapping rasters
We will use QGIS's Raster calculator. Firstly we need to create an empty georeferenced image with needed extent. Go to Raster -> Miscellaneous -> Build Virtual Raster. Here choose rasters A and B to become a virtual raster (VR) and do not define noData values. Add VR to the TOC.

Now we will create a templates for expanded A and B. In Raster calculator (Raster -> Raster calculator) insert expression: VR@1 * 0, highlight VR raster (by clicking on it) in Raster bands and click on Current layer extent button. Execute expression and save its result - it will be a template for raster A. Copy it to create template for raster B.

Now lets change extent for A. In Raster calculator enter expression:  A@1 this will simply leave A values unchanged. Then click on VR raster in Raster bands to highlight it and click a on Current layer extent button - this will set desired extent of the output raster. Save result as a template for A (overwrite existing file). Do the same for B.

Now you have both A and B with the same extent and can freely perform calculations over them.