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.

Download utility.

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.

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.

Rgdal Package: R crash during "writeGDAL" solved

Lost entire day fighting this ridicules crash during export of raster created in R to GeoTiff (other "rgdal" functions worked fine). This option worked well several months ago when I last time used it, so I wondered what happened. I searched for the only error message I was able to catch: "Failure during IO", but wasn't able to find anything useful. 

After some time I started to suspect GDAL itself in abnormal behaviour since its latest updates caused some troubles already. I was wrong but the idea itself turned to be useful - I recalled that I had about a dozen of GDAL updates since "rgdal" package was installed. Despite the package was "up to date" I've reinstalled it. It solved the problem)))

Sverdlovskaya Region, Russia: a Field Survey in Autumn. Part 2: Scars on the Earth

Ok, it's time to finish the story about land monitoring in Sverdlovskaya region. In this post I would like to demonstrate some of the most unpleasant types of the land use.

Lets begin with illegal dumping. This dump (note that there is the smoke from waste burning down) is located right next to the potato field (mmm... seems these  potato are tasty). The ground was intentionally excavated here for dumping waste. Obviously this dump is exploited by the agricultural firm - owner of this land, but who cares...

Panorama of freshly burnt illegal dump

The next stop is peat cutting. A huge biotops are destroyed for no good reason (I can't agree that use of peat as an energy source is a good one). At the picture below you can see a peat cutting with the area of 1402 ha. There are dozen of them in the study area...

Peat Cutting (RapidEye, natural colours)

But the most ugly scars on the Earth surface are left from mining works. There is Asbestos town in Sverdlovskaya region. It was names after asbestos that is mined  there. The quarry has an area of 1470 ha and its depth is over 400 meters. Its slag-heaps covers another 2500 ha... The irony is that this quarry gives a job for this town and killing it. You see, if you wand to dig dipper you have to make quarry wider accordingly. Current depth is 450 m and in projects it is over 900 m, but the quarry is already next to the living buildings. So quarry is going to consume the town... By the way, the local cemetery was already consumed. Guess what happened to human remains? Well, it is Russia, so they were dumped into the nearest slug-heap.

Here is the panorama of the quarry. You may try to locate BelAZ trucks down there ;-)

Asbestos quarry

Here is the part of the biggest slag-heap:

A slag-heap

That's how it looks from space:

Asbestos town area (imagery - RapidEye, NIR-G-B pseudo-colour composition)

And in the end I will show you the very basic schema of disturbed land in the study area (no settlements or roads included). Terrifying isn't it? 

Basic schema of disturbed land 

Conference Announcement: Open GIS!

Open GIS! is a conference for users and developers of open-source GIS. It will take place on November 17-18, 2012 in Moscow, Russia. The really cool thing is that participation in this crowd-source conference is free of charge!!! Of course modest donations will be appreciated and a sponsorship as well, but it is not mandatory - we want this event to be available for everyone who are interested in open geospatial software and data. Download leaflet!

A call for papers will be opened in August but already now you may propose topics of your talk or master-class to the committee.

Please, spread the word!

P.S. If you are (a rear) person who interested in illegal dumping research (from geospatial point of view) - come to this conference to hear my talk on it!

How To Get Rid of Needless SRTM Files Fast

It is quite annoying to download the whole 1Gb archive of SRTM30 DEMs when you need only couple of dozens of them. In case of Leningrad region for example I had to download 4 archives and face necessity of finding needed DEMs amongst several thousands(!!!) of files and deletion the rest... Too tedious and too long task to do it manually. Lets optimise it!

We will use QGIS and Linux (I prefer openSUSE) console (and a text editor of your choice). The algorithm is applicable to other similar issues too.

Firstly we will get the list of names of files that we need to save (for our actual project). Collect all DEMs from archives in a single folder. And use Image Boundary plugin (or Raster tile index from "Raster" menu) to create a polygon shp-file with borders of all the DEMs from the folder and corresponding file name and path.

Open this layer in your canvas and select polygons that intersects your area of interest:

Now you can either to "save selection as" CSV from the layer menu or "copy selected rows to clipboard" from attribute table and paste it to the empty text-file. In LibreOffice Calc (or use other text editor for this) delete all columns exept that which contains file names and add to the list name of the text-file itself (lets name it FILELIST). NOTE: this file should contain ONLY NAMES without paths. Place this file to the directory with DEMs.

We are ready for the final step - deletion of a thousands of DEMs that we are not going to use. A recipe I found here: using console go to the folder with DEMs and execute following command:

ls -1 | sort FILELIST FILELIST - | uniq -u | xargs rm

Where FILELIST is our text-file containing names of files that we need to KEEP. This command will delete all files from the folder except files which names are listed in FILELIST.

About On-line Education

Recently I've accomplished 2 of the on-line courses provided by Coursera. 

A "Model Thinking" course provided a nice overview of the existing approaches to model creation. Despite oversimplified maths (that I tried to compensate by using R-scripting where it was applicable) I liked this course for covering a wide range of real world implementation of different types of models.

A "Game Theory" course was quite nice, but had a lack of the real world implementation examples (there were a couple of them). And it was quite disappointing because that was what I hoped to learn about. I've already studied Game Theory in my University 7 or 8 years ago and since then was curious in seeing the theory in action. Unfortunately all the examples I've encountered were either oversimplified or far too complicated. The good thing was that the statement of accomplishment of Game Theory course was much fancier then that of Model Thinking.